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PREFACE 


This text is to a large extent a result of teaching two courses in molecular modeling and 
computational chemistry at the Norwegian University of Science and Technology (NTNU) 
in Trondheim. An introductory course in Molecular Modelinghas been given annually since 
2002 on the M.Sc. level based on the book by Leach (Leach 2001). This course gives an 
introduction to and an overview of the topic, including the basic elements in computational 
quantum chemistry, force fields and molecular simulations, as well as some more specialized 
topics as free-energy calculations and solvation models. A biannual course on the Ph.D. 
level, Advanced Molecular Modeling, has been given since 2004 based on own lecture notes 
and review papers. These notes have previously been used in a course organized by Prof. 
Kurt V. Mikkelsen at Aarhus University (1995) and annually at the University of Copenhagen 
(1997-2002). For this course, two sets of lecture notes, Intermolecular Interactions and 
Simulations of Liquids, were developed, where the notes on Intermolecular Interactions are 
based on an introductory chapter in my Ph.D. thesis (Astrand 1994). The notes have also 
been used in a course on Intermolecular Interactions at the University of Tromso in 2002, 
and at a summer school in Molecular Dynamics and Chemical Kinetics: Exploitation of Solar 
Energy at the University of Copenhagen annually since 2013. This text is therefore the result 
of lecture notes gathered and updated continuously over the years. 

There are many excellent books in the field of computational and theoretical chemistry, but 
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they are often specialized in one or a few of the topics important for a general course in 
molecular modeling. So it is rather the lack of a book with the, according to me, desired 
composition (table of contents) for a general text on molecular modeling than the lack of 
good texts on each of the topics that lead to that eventually this project was initiated. 

The clear separation between content and style in KTpX (Lamport 1994) makes it pivotal in 
organizing and developing a complex document. In addition, the PGF and TikZ graphics 
systems for TpX have been used extensively to construct the graphics leading to that all 
figures in the document are included as in-line LTjiX code. Consequently, the graphics 
appear in a consistent way and can be easily updated as the document is developed. All 
references with a doi are clickable in the reference lists as a result of using BibTex together 
with the doi package. Also text marked with brown (with one exception) are clickable with a 
link to an external web-page. Developing a complex LTpX document has many similarities to 
software development. Since the repository only consists of text files (including the figures 
when PGF/TikZ is used), it is therefore natural to use a version control system and for this 
project git (Chacon and Straub 2014) is used. 

There is a multitude of software available to do the actual calculations using the methods 
discussed in this text. If possible, I have so far chosen to use software that is generally 
available in the Ubuntu Linux-based system. Avogadro is used as a molecule editor (Hanwell 
et al. 2012) to generate input files for the quantum chemical calculations, and for the 
quantum chemical calculations NWChem (Valiev et al. 2010) has been used. 

There are of course many persons that have contributed indirectly to this text. I am in 
particular grateful to my Ph.D. thesis adviser Prof. Gunnar Karlstrom (Lund University) and 
to my postdoc adviser Prof. Kurt V. Mikkelsen (University of Copenhagen). Since the notes 
have been used extensively in courses over the years, I am also grateful to all the students that 
have commented on different parts of the original notes or in other ways given feedback. 

I also would like to thank Bookboon for publishing this text, and in particular I would like to 
thank Karin Hamilton Jakobsen at Bookboon for the encouragement to actually convert a set 
of separate notes into one coherent document. I support the idea of Bookboon to distribute 
free ebooks for students. 


First edition 

The 1st edition is by no means a complete book on molecular modeling, it is rather a 
compendium containing some of the chapters relevant for a general course in molecular 
modeling and computational chemistry. Apart from a brief introduction, this edition 
consists of two chapters on computational quantum chemistry and force fields, respectively. 
Since this text is published as an e-book only, it is, as for a software, possible to publish 
corrections and additions frequently. The goal is to publish a new edition annually as long 
as I use the text myself in teaching, so comments on the content are most welcome. 


POA 
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1.1 What is molecular modeling? 


With molecular models we mean models where a molecule is constituted by atoms 
connected by chemical (covalent) bonds (see figure 1.1 for two examples). Molecular models 
are normally based on a mechanistic model for describing the structure of molecules and 
other molecular properties as well as the interactions between molecules. Using the here 
quite ill-defined term atom rather than nucleus for the constituents of a molecule indicates 
that the constituent of interest consists of both a nucleus and core electrons (whereas the 
valence electrons mainly form the covalent bonds between the atoms), and we also therefore 
refer to the models discussed in this text as atomistic models. The letters in the figure denote 
the position of the atoms: H for hydrogen, C for carbon, N for nitrogen, O for oxygen, etc., 
and the solid lines denote covalent bonds, i.e. electron pairs shared by the two connected 
atoms. Another common way to depict the structure of molecules is with ball-and-stick 
models, where two molecules are shown in figure 1.2. Here a colour code is used for each 
element: black for C, red for O, blue for N, white for H, etc. 

The terms molecular modeling, computational chemistry and theoretical chemistry are often 
used interchangeably and the distinctions, if ever been meaningful, have more or less 
lost their meaning. A text on computational chemistry should in my opinion include the 
aspects of how to solve the problem on a modern computer system including the choice 
of algorithms, parallelization on large-scale clusters and optimization on gpu-accelerators. 
The subtitle of this text Concepts in Computational Chemistry indicates that we rather focus 
on what is needed from a user perspective to understand the methods in computational 
chemistry rather than the implementation of the methods. The so far never-ending rapid 
development of computer technology has evidently lead to a revolution in chemistry, and 
computational chemistry has become an accompanying analysis technique in line with 
many common experimental characterization techniques. So the role of high-performance 
computers is indisputable, but in this text we assume that we have the required computer 



Figure 1.1: Chemical structure for urea (left) and 
phenol (right). In phenol, the hydrogens on the 
phenyl ring are suppressed. 


Figure 1.2: Ball-and-stick representation of water 
(left) and formamide (right). 
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resources in terms of both hardware and software at hand. Theoretical chemistry, although 
by many often used as a synonym for quantum chemistry, is the widest term including 
theoretical development (i.e. new equations), computational chemistry (i.e. how to solve 
the equations on a computer), and how to apply the methods in different applications which 
requires a detailed knowledge of each particular application area. Theoretical chemistry is 
also a wider concept than molecular modeling and includes also thermodynamics, chemical 
kinetics and molecular informatics. 

Computer modeling is in many cases an attractive alternative to experiments with the 
requirement that the accuracy of the modeled property rivals that of the experiment. 
Computer-based simulations are in most cases less expensive and less time-consuming than 
to carry out the corresponding experiments, and in addition simulations can more easily 
be performed at extreme conditions, e.g. at very high pressures or temperatures, or with 
hazardous components, e.g. with explosives or poisonous molecules, than experiments. 
In simulations, it is also easier to investigate different contributions, e.g. from a non-zero 
temperature or from a solvent, to a property since in calculations we often add up various 
contributions which thus can be analyzed individually whereas we in an experiment often 
only get a single number as the result. Similarly, it is often possible to partition the computed 
property into various terms or in other ways to analyze the computed result in terms of 
properties that cannot be measured experimentally. As an example, the electrostatics of a 
molecule is commonly analyzed in terms of partial atomic charges, a property that cannot 
be measured experimentally. 

When presenting a method in computational chemistry, we are thus interested in three 
different things: the theory behind the method describing which properties that can be 
computed (at least in principle), the accuracy of the method, and finally, how the results can 
be analyzed to provide further insights about the studied system. To stop after the second 
step is a pity, then an accurate calculations is not more valuable than an accurate experiment 
since it only provides "a single number". 

1.2 Brief summary 

The most fundamental way to describe a molecular system theoretically is with quantum 
mechanics. In molecular quantum mechanics (quantum chemistry), we normally approx¬ 
imate both nuclei and electrons as point particles, i.e. each particle has a mass, an electric 
charge and possibly also a spin. Nevertheless, the molecular problem in quantum mechanics 
is complicated and only the hydrogen atom (one nucleus and one electron) in a clamped- 
nucleus approach (the nucleus is kept in a fixed position in space and has no kinetic 
energy) has been solved analytically. The goal is to solve the Schrodinger equation for 
molecular systems, but for many-electron atoms and all molecules this can only be achieved 
by approximate models solved numerically. A major part of computational chemistry is 
therefore devoted to approximate methods for calculating molecular energies and properties 
including very accurate molecular-orbital methods to include electron correlation, methods 
based on density-functional theory, and phenomenologically based force-field methods. If 
the wavefunction of a molecule is known, however, all information about the molecule can 
be extracted from its wavefunction. An introduction to computational methods in quantum 
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chemistry is given in chapter 2 . 

Although we briefly introduce classical electrostatics in chapter 2 , we are often interested 
in the interaction between a molecular system and an electromagnetic field. The obvious 
example is spectroscopy where the molecular system is perturbed by an applied electromag¬ 
netic field and the response is measured. Another example is organic solar cells (e.g. Gratzel 
cellswhere the light of the sun is absorbed by a photosensitizer and the excited electron is 
separated from the hole leading to an electrical current. Also long-range intermolecular 
interactions can be described in terms of electrostatics using the same basic concepts, i.e. 
a molecule is interacting with the electrostatic potential, electric field, etc. arising from the 
charge distribution of the surrounding molecules. The proper starting point of this branch 
of molecular modeling are Maxwell’s equations. 

At least historically, quantum chemical calculations are computationally too expensive to 
be used for very large systems (thousands of atoms) or in molecular dynamics simulations 
where the interatomic forces have to be computed repeatedly (perhaps millions of times). 
This gap is filled by force fields, i.e. simple and approximate models for the molecular energy 
as well as intermolecular interactions that is feasible for large-scale molecular dynamics 
simulations. Force fields are here first introduced phenomenologically and subsequently 
in a more systematic way by deriving each term in a force field from quantum chemistry. An 
introduction to force fields is given in chapter 3. 

A quantum chemical calculation gives in principle information about the properties of a 
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Microscopic world: 

Quantum mechanics 
jfey/ = Ey/ 

Classical mechanics 
F = ma 



Statistical mechanics 
S = Icb In W 



Macroscopic world: 

Thermodynamics 
U = q+w 

Kinetics 

^ = kflA]-k r lB] 



Figure 1.3: Statistical thermodynamics provides the connection between the microscopic and the macroscopic 
world. 


molecular system at the temperature 0 K. To obtain properties for a liquid or solution for 
example at room temperature and ambient pressure, we need to employ statistical thermo¬ 
dynamics. Statistical thermodynamics (statistical mechanics is used as an equivalent term 
in this text) is the theory that provides the link between the microscopic world described 
by quantum mechanics (and sometimes classical mechanics) and the macroscopic world 
described by thermodynamics and chemical kinetics (see figure 1.3). A key component of 
statistical thermodynamics is the partition function and all thermodynamics properties of 
a system can be provided from the partition function provided that it is known. This is 
not the case for a realistic system as a molecular liquid, so the problem of calculating the 
properties of a liquid is instead turned into a sampling scheme where liquid configurations 
are sampled from the correct distribution (e.g. at a given temperature, pressure and density) 
using molecular dynamics or Monte Carlo simulations. Another major part of computational 
chemistry is therefore devoted to simulations of molecular liquids. 

A chemical event, e.g. a chemical reaction or the absorption of a photon, is in most cases 
local in space where the actual event involves perhaps 5-10 atoms whereas the total system 
may consist of thousands of atoms as well as fast in time, often on the femtosecond scale. 
Consequently, multiscale and multiphysics methods in theoretical chemistry have been 
developed over the years. 

In molecular informatics, which may be subdivided into chemoinformatics and bioin¬ 
formatics, molecular properties are related to how the system functions (photovoltaic cell, 
electrochemical battery, drug, etc.) by statistics without an underlying mechanistic model 
and is therefore a separate branch of theoretical and computational chemistry. 

Many of the grand challenges in chemistry today are strongly connected to severe problems 
for our society, as for example a sustainable production of energy and electricity, clean water 
and food production, the environment, and nano-scale devices for the next generation of 
information technology. In all these cases, modeling on the atomistic scale have already 
provided or can give substantial contributions. 
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Quantum mechanics, together with statistical mechanics, is the foundation of theoretical 
chemistry and molecular modeling. Quantum mechanics applied on molecular systems, 
quantum chemistry, provides the chemical model for describing chemical bonds and 
reactions, intermolecular interactions and molecular properties. Quantum chemistry also 
provides the foundation for many less sophisticated (coarse-grained) models for describing 
molecules such as force fields. 

Previous knowledge in quantum mechanics is expected in line with an undergraduate course 
in physical chemistry (see “Recommended Literature” at the end of the chapter), and many 
of the sections in this chapter are regarded as repetition which is also reflected in the form of 
the presentation. The goal of the first sections are to provide the fundamentals of quantum 
mechanics and quantum chemistry needed in the forthcoming sections and chapters. For 
a more complete presentation of quantum chemistry, specialized texts on the subject are 
recommended at the end of the chapter. 

2.1 The Schrodinger equation 

The Schrodinger equation (Schrodinger 1926) is here presented as a hypothesis that has 
proven to be incredibly useful, and we do not aim at giving a motivation for its existence 
or the way it looks like. Knowing the solution to the Schrodinger equation provides all the 
necessary information of a microscopic system at the temperature 0 K. The time-dependent 
Schrodinger equation is given as 


a^(fi.„ JV ,f) 

dt 


t) = \h 


( 2 . 1 . 1 ) 


where H=hl2n and h is Planck’s constant, t is the time, 'Y is the wavefunction where is 
a short-hand notation for the position vectors of N particles, iy, iy ,..., Tn- The Hamiltonian, 
is the energy operator and is divided into a kinetic energy operator, ST, and a potential 
energy operator, Y, as 


je = 3- + y. 


( 2 . 1 . 2 ) 


The kinetic energy operator, ST , is the sum of the kinetic energy operator of all particles in 
the system, 



(2.1.3) 


where m,- is the particle mass of particle i, N is the number of particles, and V? is the 
Laplacian of particle i given in Cartesian coordinates as 



(2.1.4) 
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where f/ = (x/,y/,z/) is the position vector in Cartesian coordinates. The potential energy 
operator, T(ri...jv), is unique for each type of system, and in chemistry we are mainly 
interested in the Hamiltonian for molecular systems which is discussed in section 2.2. 

If the Hamiltonian is a function of only the spatial variables, and not of the time, 

separation of variables is used to simplify the Schrodinger equation. The wavefunction 
is thus written as the product of a spatial wavefunction, y/{fi ,.jv) and a time-dependent 
function, 0(f), 

v f / (Ti...iV> t) = i/dri...jv)0(f) - (2.1.5) 

which is plugged into the Schrodinger equation in eq. (2.1.1) leading to 

1 * ._ld0 

—je\fr = \ h -. ( 2 . 1 . 6 ) 

y/ 0 dt 

Since the left-hand side is a function of only and the right-hand side is a function of 
only f, both sides have to be equal to a constant, identified as the energy E. This leads to the 
time-independent Schrodinger equation for the spatial part, 

jfry/ = Ey/, (2.1.7) 

and to a trivial solution for the time-dependent part, 

-i Et 

0(f) = Ce~ , (2.1.8) 

where C is a constant. For most Hamiltonians (but not all), the solution to the time- 
independent Schrodinger equation in eq. (2.1.7) is quantized, 

= Entyn > (2.1.9) 

which is interpreted as that a quantum particle is in a state n with discrete energy levels 
at E n . Eq. (2.1.9) is an eigenvalue problem, where E n are the eigenvalues and y/ n are the 
eigenfunctions of the operator The time-independent Schrodinger equation 1 is solvable 
analytically only for a few model systems, where some of them are discussed in appendix 2.A. 


2.2 The molecular Hamiltonian 


In quantum chemistry, a molecule is represented by n electrons and N nuclei, where an 
electron has a charge —e and a mass m e , and a nucleus I has a charge Z/e and a mass mj. 
Both the nuclei and the electrons are regarded as point charges, i.e. they have no extension 
in space. The kinetic energy operator, 3~, for a molecule is given by a trivial extension of 
eq. (2.1.3) as a sum of the kinetic energy for all nuclei and electrons, 


N 


-h 2 


^ = Et^ v / + E 


-H 2 


/=! 2 mi £i 2 m e 


-Vf 


nuclei 


electrons 


( 2 . 2 . 1 ) 


the remaining part of the text, the time-independent Schrodinger equation in eq. (2.1.9) is referred to 
as the Schrodinger equation. Also, the eigenfunctions in eq. (2.1.9), y/ n , are referred to as the wavefunction. If 
the time dependence is included, we will explicitly refer to the time-dependent Schrodinger equation and the 
time-dependent wavefunction, respectively. 
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The potential energy operator, Y, for a molecule is given by the Coulomb interaction 
between point charges for all the nuclei and electrons, 


f= f | " ” -Z 7 e 2 

1 = 1 , 471£oRij i=1I=1 47t£oru 

J=I+1 v -V-' 

v -v- ' electron-nucleus 

nucleus-nucleus 


n 

+ E 


e * 


f~l, ^ n£ O r ij 
j=i +1 


electron-electron 


( 2 . 2 . 2 ) 


which is thus divided into nucleus-nucleus, electron-nucleus and electron-electron interac¬ 
tions. Here Z/e is the charge of nucleus I so that Z/ is the atomic number of the nucleus, 
e.g. Zj = 1 for hydrogen and Zj = 6 for carbon, and e is the elementary charge so that the 

charge of the electron is -e. We normally use capital letter subscripts, I,J,K .to denote 

nuclei and small letter subscripts, i, j, k,..., to denote electrons. If the distance involves only 
nuclei, it is denoted by R, whereas r is used if the distance involves at least one electron. 
In quantum chemistry, it is common to use atomic units (instead of SI units), where some 
constants are set equal to ±1, see table 2.1. In atomic units, the molecular Hamiltonian for 
the kinetic energy operator becomes 


1 iV I in 

z /=1 mi z i=1 


nuclei 


electrons 


(2.2.3) 
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charge of electron: e = -1 

length: 

lbohr = 0.529177 A 

mass of electron: m e = 1 

energy: 

1 hartree = 2625.4996 kj/mol 

H = h!2n = 1 


1 hartree = 27.2113845 eV 

4^ 

cn 

o 

II 

I — 1 




Table 2.1: Atomic units. Constants set to ±1 (to the left) and some common unit conversions (to the right). 


and for the potential energy operator we get 


II 

ZiZj 

R IJ 

f £ -Z/ 

+ LL + 

r tI 

n 

I 

/= i, 

/=7+l 


S -v-' 

j=i +1 


1 

r a 


(2.2.4) 


nucleus-nucleus 


electron-nucleus 


electron-electron 


We can thus write the molecular Hamiltonian as 

je mo[ = 3- n +3- e +f nn +y en + y ee , ( 2 . 2 . 5 ) 

i.e. the kinetic energy operators of the nuclei and electrons, respectively, as well as the 
nucleus-nucleus, nucleus-electron and electron-electron Coulomb interaction operators. 


2.3 Some basic properties of the wavefunction 


We will essentially only list some of the basic properties of the wavefunction needed in 
the forthcoming sections. For a more systematic introduction to basic molecular quantum 
mechanics, see e.g. (Atkins and Friedman 2010). 

According to Born’s interpretation of the wavefunction (Born 1926), y/*y/ i dT is interpreted 
as the probability for a particle in state i to be in a volume element dr. Here, denotes 
the complex conjugate of t/q, i.e. the wavefunction may be complex including both a real 
and an imaginary part, however the probability y/* y/ i dr has by construction only a real part 
which is a requirement for an observable. Assuming that a probability for state i, pdf), is 
normalized, i.e. the probability to find a particle anywhere in space is 1, we have 

P/dr = j y/*y/ i dr = l, (2.3.1) 

all space all space 


and we refer to this condition as if the wavefunction is normalized. Here dr is the volume 
element, which in Cartesian coordinates for a single particle is dr = dx d y dz and in spherical 
polar coordinates it is dr = r 2 sin(0) dr d6 dq>, respectively. In this text, the integration limits 
are dropped if we integrate over all space so that 



all space 


(2.3.2) 


If we return to the time-dependent wavefunction in eqs. (2.1.5) and (2.1.8) putting C to 1 in 
eq. (2.1.8), 

-i Et 

x ¥(t)=y/e^r 9 (2.3.3) 
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we see that if y/ is normalized, also Tfir) is normalized. If we in general have 


'F’OTtt) =y/*y/, 


(2.3.4) 


we refer to the state as being stationary. For an operator Cl, the expectation value, (Cl)i, is 
defined as 


<a>,- = 


/y? fly .dr 
J>*^dT 


(2.3.5) 


for a system in state i. For a normalized wavefunction, it becomes 


(n)i = 



(2.3.6) 


If y/i is an eigenfunction of n, 

f y/*Cly/. dr F2,- fw* yr . dr 

( fvfadr ~ fr t V t dr " ” 


(2.3.7) 


i.e. the expectation value is equal to the eigenvalue, QThus we can write the energy, E/, as 
an expectation value of the Hamiltonian as 


Ei = 



(2.3.8) 


for a normalized wavefunction. To simplify the notation, the Dirac bra-c-ket notation is 
introduced, 

(y/i\Cl\y/j) = (i\Cl\j) = Jy/*Cly/j dr, (2.3.9) 

where (i/q| or (i\ is the bra of state i and | y/j) or | j) is the ket of state j. The molecular 
Hamiltonian is hermitian, i.e. it fulfils 



(2.3.10) 


which leads to that its eigenvalues are real and that the eigenfunctions are orthogonal 
(see exercise 2.1 to show this). Since energies, or to be more precise energy differences, 
are measurable quantities and therefore the energies have to be real, it is a requirement 
that the molecular Hamiltonian is hermitian. For orthonormal states (i.e. orthogonal and 
normalized), 

J y/*y/j dr = {y/A y/j) = (i\j) = 6ij (2.3.11) 

where 8j j is the Kroenecker delta function (1 if i = j; 0 if i ^ j)- 


2.4 The Born-Oppenheimer approximation 

In the Born-Oppenheimer approximation, the molecular wavefunction, y/{R\...N>r\... n ), is 
approximated as the product of an electronic, y/ e \ and a nuclear, y/ nuc , wavefunction, 

~ yr d <? 1 ... n ;R 1 ...N) i^ nuc (^i...iv). (2.4.1) 
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Figure 2.1: Sketch of a potential energy surface, V{R), for a regular diatomic molecule. The minimum of V{R) 
corresponds to the equilibrium bond length of the molecule. 


where i//' el is a function of the electronic coordinates, and depends parametrically on 
the nuclear coordinates, Ri...n- 2 It means that we solve the Schrodinger equation for y/ ei for a 
given molecular geometry, the clamped-nucleus approach. The corresponding Hamiltonian, 
-i^ el , is given as 


je 


■el 


i n n N 7 n -i N 

= E E 


ZiZj 


i -1 


i=l/=l r U 


i =1 r ij 
j=i +1 


/=1 R IJ 
/=/+! 


— 3~ P + 7 gfi + Tgg + T^nn > 


(2.4.2) 


i.e. the kinetic term for the nuclei is ignored in the clamped-nucleus approach as compared 
to the molecular Hamiltonian in eq. (2.2.5). The last term on the right-hand side in eq. (2.4.2), 
the nucleus-nucleus potential energy, becomes a “constant” contribution to the molecular 
energy since the nuclear positions are regarded as parameters and not as variables in the 
electronic wavefunction i// el . The Schrodinger equation for the electronic state i becomes 

= £ ft. (2.4.3) 

If eq. (2.4.3) is solved repeatedly for different molecular geometries, a potential energy 
surface is obtained for each state i, ef{Ri...N)- Normally, we refer to the ground state energy 
surface, £o 1 (-Ri...jv)> as the potential energy surface, V(Ri ^j), 

VtRi...N) = £%\Ri...n) , (2.4.4) 

where a typical potential energy surface for a diatomic molecule is depicted in figure 2.1. 
The zero-level of the energy scale is normally shifted for FCRi..jv) as compared to cf(R]...N) 
so that V{Ri_n) approaches zero for an infinite separation of two fragments, whereas 
£q 1 (-/?i...]v) approaches zero for an infinite separation of all nuclei and electrons. We define 
the Hamiltonian for the nuclei as 

N i 

& nuc = - E ^- V / + V ( R 1 ...n) > (2.4.5) 

j=i 2m/ 

2 Note the distinction between /(x,y) and /(x;y). In the first case, x and y are both variables, but in the 
second case, y is a parameter, i.e. it has a single, constant value. 
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and the corresponding Schrodinger equation becomes 

^nucynuc = £ nuc^nuc _ (2. 4 .6) 

By applying the Born-Oppenheimer approximation, we have thus separated quantum 
chemistry into two problems: the electronic problem in eq. (2.4.3) solved for a given 
molecular geometry, and a nuclear problem in eq. (2.4.6) where the potential energy surface 
is obtained by solving the electronic Schrodinger equation for a set of molecular geometries. 
In this chapter, we focus entirely on the electronic structure of molecules by discussing 
methods for solving eq. (2.1.8). In a forthcoming chapter on molecular structure and 
vibrational motion, we will discuss how to solve eq. (2.4.6). 


2.5 Atomic orbitals 

2.5.1 One-electron atom 

The starting point for solving the electronic Schrodinger equation is the one-electron 
atom. The position of the nucleus is regarded as fixed by adopting the Born-Oppenheimer 
approximation in section 2.4. The Hamiltonian, <^? el , for an electron interacting with a 
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/ 

0 1 

2 

3 

4 

5 

name 

5 p 

d 

/ 

g 

h 


(a) Symbols used for the quantum 
number l 


n 

/ 

mi 

name 

1 

0 

0 

Is 

2 

0 

0 

2 s 

2 

1 

0,±1 

2 P 


(b) Solutions allowed for 
n - 1,2 


Table 2.2: Notation for atomic orbitals 


nucleus becomes in atomic units, 

1 Z 

,i£ el = —V 2 --, (2.5.1) 

2 r 

where Z is the charge of the nucleus and r is the distance between the electron and 
the nucleus. The Hamiltonian thus consists of a kinetic energy operator for the electron 
and the Coulomb interaction between the electron and the nucleus. The solutions to the 
Schrodinger equation is in spherical polar coordinates given as (see appendix 2.A.5) 

nlmi } (p) — Rnl (t") ZifYii [0 , (p) , (2.5.2) 

where f?„/(r) is a radial function and Yi mi (6,(p) is a spherical harmonics. The solution 
depends on three quantum numbers, the principal quantum number n and two angular 
quantum numbers l and ra/, which are restricted to the following integer values: 

n = 1,2,3,... 

/ = 0,1,2,..., n- 1 

m i = 0 , + 1 , + 2 +1 

The naming convention for atomic orbitals is given in table 2.2, giving the notation of Is, 
2s, 2 p, etc. orbitals. We thus have three 2p-functions, normally denoted 2 p x , 2p y and 2 p z . 
For n = 3, we get 3s, 3 p (with the components 3 p x , 3 p y and 3 p z ), and 3d functions. The 
^-functions have five components, normally labeled d xy , d xz , d yz , d z 2 , and d x 2 _y>. It is also 
noted that the atomic orbitals form an orthonormal set of functions, 

J Wnlmi dT — . (2.5.3) 

Each electron has a spin, with a spin quantum number, m s , 

m s = ± (2.5.4) 

Each spatial atomic orbital in eq. (2.5.2), may thus accomodate two electrons without 
violating the Pauli principle for fermions (each electron has a unique set of n, l, mi, and 
m s ). The electron configuration for an atom is given according to the Aufbau principle as for 
example 

He Is 2 
Ne ls 2 2s 2 2p 6 

Cl ls 2 2s 2 2p 6 3s 2 2p 5 orNe3s 2 2p 5 


Download free eBooks at bookboon.com 


12 








ATOMISTIC MODELS 


MOLECULAR QUANTUM MECHANICS 


2.5.2 Two-electron atom 


If we next regard the Hamiltonian of the two-electron atom, still within the clamped- 
nucleues approach, 





“V— 

A 


z 1 , z 1 

- v 2 i - —+— > 

'•/ 2 1 '•./• r ij 

-V-' 

Aj 


(2.5.5) 


where we have two electrons i and j interacting with a single nucleus. The Hamiltonian 
consists of five terms: a kinetic energy for each electron, the Coulomb interaction between 
each electron and the nucleus, and the repulsive Coulomb interaction between the two 
electrons. If we as a first approximation ignore the electron repulsion, the Hamiltonian 
becomes 

,76^ = ,Ai + , (2.5.6) 


i.e. two one-electron Hamiltonians of the type in eq. (2.5.1) each with a solution given by 
eq. (2.5.2). Denoting a one-electron wavefunction, an atomic orbital, with (piifi), we may 
anticipate that the Schrodinger equation becomes using variable separation 


[A + Wj) (pi ( 7i)<pj ifj) = (e, + £/) ^C fi)(Pj ( fj ) , (2.5.7) 

where c* is an orbital energy. Since electrons are fermions, they are, however, indistin¬ 
guishable and the wavefunction has to be anti-symmetric with respect to the exchange of 
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two electrons. If both electrons are put into the Is orbital, an occupation denoted Is 2 , the 
wavefunction may as a first attempt be written as 

= ls(z) ls(j), (2.5.8) 

where the notation f/ = i and Tj = j is adopted. In eq. (2.5.8), the electrons are 
indistinguishable but the wavefunction is not anti-symmetric. It is the total wavefunction, 
however, including also the spin in addition to the spatial coordinates, that has to be anti¬ 
symmetric. Denoting a spin-orbital, Xi(^i> h)> 


i’^i) — (piifi) & ii$i) > 


(2.5.9) 


where a, can be either a for m s = \ or for m s = However, the wavefunction, 

ys(i,j) = 1 s(i)a(i) x 1 s(y) /3(y), (2.5.10) 


is not anti-symmetric (nor are the electrons indistinguishable). Constructing the following 
linear combinations of the spin-functions, 


1 

7 ^ 

i 


(a(i)/3(j') + a(y')j8(t)) 
[a{i)P{j)-a{j)fKi)] 


symmetric, 
anti-symmetric, 


gives a symmetric and an anti-symmetric spin-function. Here it is assumed that the spin- 
functions are orthonormal, i.e. 

(a(i)\P(j)) = 6 a p8 ij , (2.5.11) 

giving the normalization factor, IIV2, above. The wavefunction, 

ys(i,j) = ls{i)ls(j) x 7— (a(/))6(y) - a{j)(5{ij) , (2.5.12) 

v2 


is thus acceptable since the wavefunction is anti-symmetric and the electrons are indistin¬ 
guishable and in line with our chemical picture that the helium atom has a Is 2 configuration 
with a spin-up and and a spin-down electron given schematically in figure 2.2a with two 
electrons in the orbital with the lowest energy. For an excited state of He, e.g. ls(l) 25(2), see 
figures 2.2b and (c), the spatial part of the wavefunction becomes 


7- (ls(z)2s(y) + ls(j)2s(i)] , 

v 2 


(2.5.13) 


where the plus sign gives a symmetric function and the minus sign gives an anti-symmetric 
function. If they are then combined with appropriate spin functions, anti-symmetric 
wavefunctions can be constructed for both the singlet state and the triplet state. 


Download free eBooks at bookboon.com 


14 


ATOMISTIC MODELS 

1 

MOLECULAR QUANTUM MECHANICS 

T 

-Hr 

h — 

1 

-i — 

(a) Ground 

(b) Singlet 
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excited 


state 
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Figure 2.2: Some electronic states for the He atom 


2.5.3 n-electron atom 


Generalizing to an n-electron atom, a Slater determinant (Heisenberg 1926, Dirac 1926, 


Slater 1929), 



*iU) 

*i(2) 

*i(«) 


*2(D 

* 2 ( 2 ) 

*2 (/*) 


Xn( 1 ) 
XnC 2 ) 

*n(n) 


(2.5.14) 


gives the correct symmetry of the wavefunction as expressed by spin-orbitals. 


2.6 Molecular orbitals 

A fundamental approximation in molecular orbital theory is to construct molecular orbitals, 
(pi, as a linear combination of m atomic orbitals, <pj, the LCAO approximation, 

m 

cpi = E CijcPj , (2.6.1) 

j =i 

where Cij is an orbital coefficient. For the hydrogen molecule, H 2 , a molecular orbital, oy 
(where a indicates that it is a cr-bond), may be written as a linear combination of two atomic 
orbitals, a Is orbital for the first hydrogen atom, ls^, and a Is orbital for the second hydrogen 
atom, 1 s B , 

Gi = C iA 15^ + C iB l SB • (2.6.2) 

Using only two atomic orbitals for the hydrogen molecule, the solutions may be found 
trivially by regarding the symmetry of the molecule. If the two nuclei are placed at (x e ,0,0) 
and (-x e ,0,0), respectively, the probability to find an electron in x and -x will be equal 
because of symmetry reasons. Assuming real wavefunctions, 

y/ 2 (x)=y/ 2 {-x) , (2.6.3) 

which has two solutions, a symmetric and an anti-symmetric solution, 

i/dx) = tyd-x) and y/(x) = -y/(-x). (2.6.4) 
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Figure 2.3: Molecular orbital diagrams of the hydrogen molecule 


Thus, 15^ and 1 sb are two identical functions (atomic orbitals), centered around x e and -x e , 
respectively. Furthermore, cai - csi or cai - - cbi■ If the atomic orbitals are orthonormal, 
see eq. (2.5.3), orthonormal molecular orbitals are given by 

1 

where a g (the subscript g denotes gerade, German for even) is a bonding cr-orbital and a u 
(the subscript u denotes ungerade, German for odd) is an antibonding cr-orbital. Following 
the Aufbau principle, the orbital with the lower energy a g will be doubly occupied (two 
electrons with opposite spin), whereas a u is unoccupied. 

In general, the molecular wavefunction may be written as a Slater determinant of molecular 
orbitals. For a molecule with n electrons, | molecular orbitals are occupied (two electrons 
with opposite spin in each molecular orbital), whereas the remaining molecular orbitals are 
unoccupied. As an example, the ground state of the hydrogen molecule, y/o, may thus be 
written as a Slater determinant of the occupied spin-orbitals, 


(15a-15 B ), (2.6.5) 


(Jg = —[Is a + Isb) and a u 

V2 


1 Zld) I2d) 
s/2 Il(2) 12 ( 2 ) 


1 


(ll (1)12(2)-U (2)12(1)) , 


( 2 . 6 . 6 ) 


where %\ ( j) = a g (j)a(j) and ^ 2 ( 7 ) = °g( j)P(j) and the wavefunction is normalized. 
Molecular orbitals may be depicted in molecular orbital (MO) diagrams, given for the 
hydrogen molecule in its ground state in figure 2.3a. If we regard the triplet state of the 
hydrogen molecule, the MO diagram is given in figure 2.3b, the wavefunction, y>\, becomes 


1 

Zi(l) 

Tad) 

1 

<T g (l)a(l) a M (l)a(l) 

n 

Ti(2) 

Ta(2) 

~^2 

cr g (2)a(2) a u (2)a(2) 


where we have used the notation XsW = cr u (j)a{j). 
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2.6.1 Energy of the hydrogen molecule 

The electronic Hamiltonian for the hydrogen molecule is 

a pi 1 ? 1 p Z A Z A Zb Zb 1 ZaZb 

je el = —v? - -v§ - — - — - — - — + — + -, 

2 2 via r 2 A n B r 2B r i2 Rab 


( 2 . 6 . 8 ) 


which we rewrite as 


H2 


ZaZb 
Rab ’ 


where the one-electron term, J6i, for electron i is given as 


(2.6.9) 


~ 1 ? Za Zb 

je t = —v? - — - —. 

2 * r iA r iB 


( 2 . 6 . 10 ) 


We introduce 

E 0 = Eoil) + E 0 {2) + EoO.,2) + , (2.6.11) 

Rab 

where the last term on the right-hand side is a classical Coulomb energy between the two 
nuclear charges. For the one-electron contribution to electron 1, we get for the ground 
state of the hydrogen molecule in eq. (2.6.6), using the expectation value for the energy in 
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eq. (2.3.8) (see exercise2.3), 

£o(D = (V'ol'&iW'o) 

= ^ (<Il(l)l^llll(l)> + <l2(l)l^lll2(D>) 

= <C7 g (l)|^i|0- g (l)). 

Analogously, 

£ 0 (2) = <cr g (2)|^ 2 |o-g(2)> . 

The two-electron term, E 0 { 1,2), becomes (again see exercise 2.3) 

£ 0 U,2) = <cr g (l)(7 g (2)| — |cr g (l)cr g (2)) , 

12 


which may be rephrased as 


£o(l,2) 


/ 


cr g (l)o- g (2)-^-o- g (l)cr g (2) dr 



( 2 . 6 . 12 ) 

(2.6.13) 


(2.6.14) 


(2.6.15) 


where we in line with Born’s interpretation in eq. (2.3.1) relate the electron density for a 
molecular orbital as p(j) = a* ( j)cr(j). Thus, £(1,2) is interpreted as a Coulomb interaction 
between two charge distributions and these Coulomb terms are normally denoted as hj . We 
therefore rewrite eq. (2.6.11) as 

ZaZp, 

E 0 = H 1 + H 2 + J 12 + -4-1 , (2.6.16) 

EaB 

where we also have adopted the notation Hj = E 0 {j ) for the one-electron terms. 

In addition, the energy of the triplet state of the hydrogen molecule, see eq. (2.6.7) for its 
wavefunction, is obtained (see exercise 2.3). The one-electron terms become analogous to 
the ground state, 


Ei(j) = ^ [(a g (j)\^j\a g (j)) + (a u U)\^j\auU))) , (2.6.17) 

but the two-electron term, £i(l,2), becomes 

£i(l,2) = ]- f<a g (l)(7 M (2)| — |a g (l)a u (2)> + <a M (l)<r g (2)| —|a M (l)<r g (2)>] 

2 V r 12 r 12 ) 

v -v-' 

Jl2 

- <<7 g (l)<r„(2)| — |a g (2)<7 M (l)> . (2.6.18) 

12 

'-V-" 

Ki2 

Again, the first term on the right-hand side is interpreted as a Coulomb interaction, Jij, 
whereas the second term is referred to as an exchange integral and is denoted as Kij. The 
energy for the triplet state of the hydrogen molecule thus becomes 


Ei - H\ + H 2 + J 12 ~ Kn + 


ZaZb 

R-ab 


(2.6.19) 
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Since an exchange integral, Kjj, has a positive sign, it means that the triplet state in eq. (2.6.7) 
has a lower energy than the corresponding excited singlet state, 



crg(l)a(l) cr m ( 1)/3(1) 
crg(2)a(2) a M (2)/3(2) 


( 2 . 6 . 20 ) 


and the interpretation is that the electrons are more delocalized when the system is in the 
triplet state. 


2.6.2 Energy of a Slater determinant 

The Slater determinant in eq. (2.5.14) is rewritten as 


( 2 . 6 . 21 ) 


where si is the antisymmetrizing operator generating the correct Slater determinant by 
operating on the direct product of spin-orbitals, si is given as 


l n_1 l 

j=—— y = — 

Vn\ p^o VrA 


n n n n n 

i-E E &if + E E E 

i=lj=i+l i=lj=i+lk=j+l 


( 2 . 6 . 22 ) 
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where 1 is the identity operator, is a permutation operator permuting the coordinates 
of two electrons i and j, 

&lf\XiWXjW) = I XjWXiU)) ■ (2.6.23) 

Analogously, gives all possible permutations of the coordinates of three electrons, 

&lf k \Xi(i)XjWXk(k)) = \Xk(i)XiWXj(Q) + \Xj(i)Xk(j)Xi(Q) , (2.6.24) 

etc. By convention, we order the orbitals in the direct products in eqs. (2.6.23) and (2.6.24) 
by the label of the electronic coordinates. It can be shown that si commutes with J€, 


sd,je 


sdje-jesd = o, 


(2.6.25) 


and that 


sisi = \fn\si. 


(2.6.26) 


Eqs. 2.6.25 and (2.6.26) are derived in exercise (2.4). We rewrite the electronic Hamiltonian 
in eq. (2.4.2) in line with eq. (2.6.9) as 


je* x = &- e + f en + f ee +f nn = 


n n n 


Lki)+Z E §(i> j) + inn > 

i =1 i—1 j=i +1 


(2.6.27) 


i.e. the one-electron term h{i) is the motion of electron i in the potential of all the nuclei 
and includes thus the 3~ e and Y ne terms and g{i,j) is the two-electron term including the 
electron-electron Coulomb repulsion term Y ee . The energy of the Slater determinant in 
eq. (2.6.21) is given as 


E 0 = <i/d^ e V> 

= (sixi ( 1 ) 12 ( 2 )... Xn(^)l^ e Wli(l)T 2 ( 2 )...I„(n)> 

= Y^.(X 1 (UX2 (2) • • -Xnin) \J^\J X 1 (1)12 (2) • • • Xn (»)> 
n -1 

= E (-l) p (XiWX2i2)...Xnm^ >el \^ ip] XiWX2i2)...Xn(n)). (2.6.28) 

p =0 

The nucleus-nucleus Coulomb operator, Y nn , only depends on the nuclear coordinates, 

(V/\inn\V/) = Vnn = Vnn > (2.6.29) 

where we have used that if/ is normalized. V nn is thus reduced to a classical Coulomb 
interaction energy as in eq. ( 2 . 6 . 11 ) for the hydrogen molecule. Using that the spin-orbitals 
form an orthonormal set, only the identity operator 1 in the expansion of si in eq. ( 2 . 6 . 22 ) 
gives a contribution to the energy of the one-electron operator h{i), e.g. 

<Ii (1)12(2)... Xnin)\ha)\xiil)X2&)...Xnin)) 

= (XiW\kl)\XiW)(X2{2)\X2{2))..AXn{n)\Xnin)) 

= <Ti(l)|h(l)l;£i(l)> = h\ . (2.6.30) 
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For the one-electron operator, all energy terms including a permutation gives zero, e.g. 

(XiWX2(2)...Xn(n)\h{l)\&^xiWX2(2)---Xn(n)) 

= (xia)m)\x2a))(x2i2)\xi(2))...(xnmxnin)) = o, (2.6.31) 

where the second term on the right-hand side with the integration over the coordinates of 
electron 2 is zero because of the orthogonality of the spin-orbitals 1 and 2. For the same 
reason, only the identity operator 1 and the two-electron permutation operator 2? a) in 
eq. (2.6.22) give contributions for the two-electron operator g{i, j ). The term for the identity 
operator for electrons 1 and 2 becomes 

<Tid)T2(2)T3(3)...T M (n)|g(l,2)lxi(l)T2(2)T3(3)...T M (n)> 

= <Ti(l)l2(2)lg(l,2)lxi(l)l2(2))<l3(3)ll 3 (3)>...<T w (n)|^(n)> 

= <Ti(l)l2(2)lg(l,2)lxi(l)l2(2))=/i2, (2.6.32) 

and is referred to as a Coulomb integral in line with eq. (2.6.15) for the hydrogen molecule. 
The second term for becomes 

<X1 (1)12 (2)13 (3)... Xn in) |g(l, 2) \0>$X1 («^2 (2)* 3 (3) • • • Xn in)) 

= <Tid)l2(2)lg(l, 2)lx 2 (l)li (2))<i 3 (3)ll3(3))...<T„(n)li„(n)> 

= <Ii(l)T 2 (2)lg(l,2)IX2(l)Ti(2)) = r 12 , (2.6.33) 

where K\ 2 is denoted an exchange integral in line with eq. (2.6.18). The combination of Slater 
determinants and orthonormal orbitals to reduce the molecular energy into a sum of one- 
and two-electron integrals is referred to as the Slater-Condon rules (Slater 1929, Condon 
1930). In eqs. (2.6.32) and (2.6.33), we put electrons 1 and 2 in orbitals 1 and 2, however 
we integrate over the electronic coordinates so the electrons could have had any labels. We 
can therefore drop the electron labels and write the Coulomb and exchange integrals as 

Ju = (XiX 2 \g\XiX 2 ) and K 12 = (XiX2\g\X2Xi'> > (2.6.34) 

respectively, but where we now have to be observant on the order of the orbitals in the 
integrals. The energy in eq. (2.6.28) may thus be written as 

Eo = y hi + ^ ^ [jij — Ki j | + V nn , ( 2 . 6 . 35 ) 

1=1 1=17=1+1 

where the minus sign for the exchange integrals arises from the i~l) p factor in eq. (2.6.28).. 
Utilizing that the self-interaction Ja is exactly cancelled by Ka, see eqs. (2.6.32) and (2.6.33), 
eq. (2.6.35) is rewritten as 


Eo = ±h i+ l -±± [jij - Kij) + V nn ■ (2.6.36) 

i =1 2 i=1 j —\ v ' 

For a closed-shell system, we put two electrons in each spatial orbital but with opposite spin. 
As an example, we take 


*i(l) = Vi(l)a(l) and * 2 (2) = <Pi(2)/3(2) (2.6.37) 
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Let us compute the exchange integral in eq. (2.6.33) for this pair of orbitals, 

Kn = (Xi(l)X2(2)|g(l,2)| Zl (2)x 2 (l)> 

= <^i(l)^2(2)|g-(l,2)|^ 1 (2)^2(l)><«(l)l^(l))<«(2)|/3(2)> = 0 (2.6.38) 

which becomes zero because of the orthogonality of the spin functions, see eq. (2.5.11). 
Therefore half of the exchange integrals in eq. (2.6.36) vanish so that 

n/2 n!2 nl2 

E 0 = 2 £ ( 2 Jij - K U - (2.6.39) 

i =1 i =1 j =1 V ' 

where the sums now run over the number of doubly occupied orbitals. 


2.7 The variational principle 


Variation theory is a method to obtain approximate solutions to the Schrodinger equation. 
We denote the exact wavefunction with i/q and an approximate trial function by i/q. The 
exact eigenenergies are denoted Ei. The energy of the trial function, Ei, is written as an 
expectation value, 


~ _ (y,-|«j£|yq) 


(2.7.1) 


which is termed the Rayleigh ratio. The variation theorem states (shown in exercise 2.5) 


E 0 > E 0 for any choice of t/q , 


(2.7.2) 


where the equal sign holds only if the trial function is equal to the exact wavefunction. 
Consequently, the energy serves as a measure of how good the trial wavefunction is, and 
we search for a trial function with the lowest possible energy. An important use case of the 
variational principle is when the trial function is expanded in a set of m functions, <p p , 

m 

{//q = Cptpp , (2.7.3) 

P =i 

where c p is a coefficient to be determined. The Rayleigh ratio in eq. (2.7.1) becomes 

m 

H CpCqHpq 
p,q= 1 

Eq = P — -, (2.7.4) 

L CpCqSpq 
p,q=l 


where we have used the notation H pq = {(f)p\jfe\<pq) and S pq = (<p p \(f) q ), respectively. 
Applying the variation theorem in eq. (2.7.2) to this trial function leads to the following 
conditions 


dEo 

dc r 


= 0 


Vr . 


(2.7.5) 
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The result is obtained as a secular equation (see exercise 2.6 ) 

TYl 

Y,c p (H pr -EoS pr ) = 0 Mr , (2.7.6) 

P=i 

which is fulfilled if its secular determinantis zero (also see exercise 2.6 ), 

|H-£oS|=0, (2.7.7) 

where H pr is a matrix element of H and S pr is a matrix element of S. This method is termed 
the Rayleigh-Ritz method. 


2.8 Perturbation theory 

2.8.1 Time-independent perturbation theory 

Rayleigh-Schrodinger perturbation theory (RSPT) is introduced (see e.g. (Hirschfelder et al. 
1964)). The purpose is to solve the eigenvalue problem for the Hamiltonian, Jt, 

Je\y k ) = + A^i) \y/ k ) = £k\Vk )» (2.8.1) 
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where the Hamiltonian is divided into two parts, and . It is assumed that the solutions 
are known for the unperturbed Hamiltonian 

( 2 . 8 . 2 ) 

whereas is regarded as a perturbation and A is an order parameter. It is thus implied that 
in some sense is small compared to and thus that the zeroth-order wavefunction, 
y/^\ is also relatively close to the exact wavefunction, y/^. Both the wavefunction, y// t -, and 
the energies, ejt, are expanded in A as 

c fc = cj™ + Ac™ + A 2 £™ + ••• and \y/ k ) = \y/ { ° ] ) + A\y/^ ] ) + X 2 \y/f) +... , (2.8.3) 

where \ty^) and c^” 1 are the nth order corrections to the wavefunction and the energy, 
respectively. We proceed by putting eq. (2.8.3) into the Schrodinger equation in eq. ( 2 . 8 . 1 ), 


(Mb + A^i) (|i/4 0) > + A|^™> + AV™) + ...) = 

= (c® + Ac™ + X 2 ef +...) [\yrf) + Ah/r™) + A V® > +...) . (2.8.4) 

This equation must hold for each order n in A” leading to for the lowest orders: 

n = 0 Wotyf) = ef\yrf) (2.8.5) 

n = 1 = c^V™) +£ fc ) l^ 0) > (2.8.6) 

n = 2 'X’o\yrf)+M’i\V% ) ) = +£^ ) \if/^ ) ) + ef ) \iff^ ) ), (2.8.7) 

i.e. the leading term in the expansion is the zeroth-order solution in eq. (2.8.2). The 
normalization of y/^™’ is chosen according to intermediate normalization as 

<*r <>=*-• (2 - s - 8) 

The energy for each order n are obtained by projecting eqs. (2.8.5)—(2.8.7) by then 

using eq. (2.8.2) and applying the condition in eq. (2.8.8), which leads to 

n = 1 c™ = (y/f |^|y/™> (2.8.9) 

n = 2 ef = <^ 0) |^il^™>. (2.8.10) 


The first-order correction to the energy, c™, is given as the expectation value of the 
perturbation operator, Jt\. For obtaining a first-order correction to the wavefunction in 
Rayleigh-Schrodinger perturbation theory, it is expanded in the spectrum of the unper¬ 
turbed wavefunction, 

OO 

Wk ) ) = 'L cl ij\v'f)> (2-8.li) 

and by projecting with (y/ 0) |, we obtain 

<«) = c™. (2.8.12) 
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The intermediate normalization condition in eq. (2.8.8) gives that cj^ = 0, which leads to 


W [ k) = £ • (2.8.13) 

jVJfc 

Furthermore, eq. (2.8.6) is rewritten as 

(4® i<> = (•*! - 4”) *?> • < 2 - 8 - 14 ) 

which is projected by (i// 0) | giving 

^ 0) - X 0) ] (y/ 0) \y/^) = |.i^i | ty[ 0) ) for all j apart from j = k. (2.8.15) 

Thus c^j in eq. (2.8.12) becomes 


AD 

"kj 


(yfmy/f) 

,40) _ AO) 

£ k £ j 


for j ^ k. 


(2.8.16) 


The first-order correction to the wavefunction in eq. (2.8.13) may thus be written as 


co <y/ 0) |.i^r> 

<> = E J k 

jfk 


AOK 


_co) _ p (0) 

k j 


l< } > • 


(2.8.17) 


Eqs. (2.8.10) and (2.8.17) are combined to give the second-order contribution to the energy, 




(0)\/,,,(0)i 


A0K 


4 2) = 


j^k 


J 


J 


-CO) _ ro) 
£ k £ j 


(2.8.18) 


Also the second-order correction to the energy, £“ , is thus given in terms of the solution 
to the zeroth-order problem in eq. (2.8.2). From eq. (2.8.18), it becomes, however, apparant 
where Rayleigh-Schrodinger perturbation theory fails. When e|, 0) « e®\ eq. (2.8.18) diverges. 


2.8.2 Time-dependent perturbation theory 

We consider a time-dependent Hamiltonian, ^{t), partitioned as 

+ Ao^i (f) (2.8.19) 

where J€\ ( t) is a time-dependent perturbation and A is an order parameter as in eq. (2.8.1). 

is the time-independent Hamiltonian for the unperturbed system for which the 
wavefunction is stationary (see eq. {2.3.4)), 

-iA 0) t 

= |i/4 0) > = e~ ia)nt |i/4 0) > (2.8.20) 
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where we introduced a> n = lh in the last step. The solution to the time-independent 
Schrodinger equation in eq. (2.1.9) is known for the Hamiltonian J&q, 

^oll/4 0) > = 4°Vn 0) > (2.8.21) 

where forms an orthonormal set (see eq. ( 2 . 3 . 11 )), 

(V'm\V'n ) )=8mn- ( 2 - 8 . 22 ) 

We also assume that before the perturbation is turned on (chosen as t = 0), we are in state p 
so that 

1^(0)) = tyf) • (2-8-23) 

The evolution in time of the wavefunction 'T(t) is given by the time-dependent Schrodinger 
equation in eq. ( 2 . 1 . 1 ), 

. (2.8.24) 

at 

We proceed by expanding the wavefunction T / (f) in the solutions of the unperturbed system, 

OO OO 

mt)) = X a n {t)\^\t)) = £ a n {t)e~ io)nt \y/^ ] ) , (2.8.25) 

n =0 n —0 
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where we used eq. ( 2 . 8 . 20 ) and a n {t) are time-dependent coefficients to be determined. We 
combine eqs. ( 2 . 8 . 25 ) and (2.8.24) leading to 


•&{t) Y a n {t)fT ia>nt \y { n ] ) =ih— Y a n it)e~ 1Wnt |i/4 0) > 


72=0 


dt ^o 


(2.8.26) 


Using eq. (2.8.19) and the dot notation (Newton’s notation) for time derivatives, 

d a n {f) 


Unit) 


d t 


(2.8.27) 


we get 


°° , . 

Y [ £ n ]a n(t)e~ 1C0nt |i/4 0) > + e~ lc ° nt |i/4 0) >) = 

n= 0 V ' 

Y [iha n {t)e-^ 1 \xffW) + efa n {t) |^ 0) >) 


72 = 0 


where two of the terms cancel each other. Projecting by 


gives 


rt(f) \ = Wm\^ Mmt 


( 0 )\ aitomn t — 


Y a nW<Wm I (■01' V [ n > e 

72=0 


= iha m {t), 


(2.8.28) 


(2.8.29) 


(2.8.30) 


where - a) n and we have used the orthonormality condition in eq. ( 2 . 8 . 22 ). We 

expand the time-dependent coefficients a n {t) in the order parameter A, 


00 

a n {t) = Y ^a^it). 

i =0 


(2.8.31) 


After substituting eq. (2.8.31) into eq. (2.8.30), the resulting equation has to hold for each 
order of A l ,i = 0,1,..., 00 . For i = 0, we get the unperturbed and time-independent 
wavefunction so a^ 0) = 8 pn from eq. ( 2 . 8 . 23 ). For i = 1, we get 

OO 

iha™{t) = Y l^iU)|i/4 0) >e iaW = <^™|^i(t)|^ 0) >e iWm P f , (2.8.32) 

72 = 0 

where we in the last step have used that aJJ” = 8 pn . We can integrate this equation, assuming 
that the perturbation is turned on at t = 0, as 


a^(f) - a 


(i) 


( 0 ) = yJ(V m \y/ p ^)e it0mpt 'dt ’, 


(2.8.33) 


where < 2 ^ (0) = 0 is given by the initial condition at t = 0 in eq. (2.8.23) and a,„ (0) = 8 mp . If 
(f) is of the form = T fit) where Y is independent of time and/(f) isafunction of 

time only, 


(f) 


= -£<Vm\y\Vp ) ) f /(f 7 ) e i( ° mpt ' dt’ . 


(2.8.34) 
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Figure 2.4: The probability versus co. It is noted that the central peak grows as t 2 and narrows as 1 It. 

If the wavefunction in eq. (2.8.25) is normalized, 

/ oo 

'P*(tmf)dT= E \ a m(t)\ 2 = 1, (2.8.35) 

m =0 

where \a m {t)\ 2 has the interpretation of being the probability, P p ^ m (t), of the system being 
in state m at time t if it was in state p at t = 0, 

Pp^ -m(f) = I • (2.8.36) 

Further developments depend on the form of /(f), where examples are Fermi’s golden rule 
in section 2.8.2.1 and the frequency-dependent polarizability in section 2.9.4. 


2.8.2.1 Fermi’s golden rule 


If it is assumed that the perturbation, (f), is independent of time apart from being turned 
on at t = 0, we have from eq. (2.8.34) that /(/) = 1 leading to (see exercise 2.13) 


flm(f) = 


(Vm\y\^p) -l 


h 


CO 


m^p. 


mp 


(2.8.37) 


Utilizing eq. (2.8.35) to first order in A, we get (see exercise 2.14) 

* 2 1 

d( 1) m - I„(l)mi2 -/HA„(0)n/h„(0)M2 Sm h ^mpt) 


l {t) = \a^{t)\‘ L = 4\(y™\y\y™)V 


h 2 coi 


m^p, 


(2.8.38) 


mp 


which is the first-order probability for the transition from state p to state m. This probability 
is shown as a function of co mp in Figure 2.4, where the most likely transitions are to states 
whose energies lie under the central peak. Since this peak is given by the first zeros of sin(x), 
the energies of the most probable states satisfy 

^ 2nh 

hco mp = \E p -E m \< —— . (2.8.39) 
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Furthermore, it is supposed that there is a continuum of states around state m, and that we 
would like to know the probability of the transition into this group of states rather than into 
a single state m. Let us assume that \(y/m is a smooth and slowly varying function 

of m and let p{E m ) be the density of states around E m . The probability of a transition from 
state p to a state around m becomes 


Ppl m W = \(Vm 



d E m 


(2.8.40) 


where d E m = h da> m and the asterisk indicates that we integrate only over a small region 
around E m . For large enough t, the central peak will include all the states around E m . If 
we also assume that p{E m ) is a slowly varying function of E m , we get 


P$l m tt) = \{y™\y\yjf)\ 2 p{E m ) J 

* 



d£m 


(2.8.41) 


For large t the area under the central peak is essentially all the area and we can extend the 
limits to +oo and use 


/ 


sin 2 (x) , 

--—d x=n 

x z 


(2.8.42) 


-oo 


The probability is written in terms of the transition rate T as 


p£l m (t) = rt 


(2.8.43) 


so that 

T = ^-\(^\f\^ ] )\ 2 p(E m ), where \E p -E m \<^ } (2.8.44) 

which is known as Fermi’s golden rule. A similar result may be obtained by noting that as t 
becomes large the probability in figure 2.4 becomes more and more peaked around E m = E p . 
Since it has the total area 2ntlh, it therefore approaches ^ 5(E p - E m ) apart from negligible 
oscillations in the wings. We thus have 

27T 

T = —\(y/^\f\y/^)\ 2 8lE p -E m ) (2.8.45) 

as an alternative form of Fermi’s golden rule. 


2.8.3 Use cases for perturbation theory 

In this text, we will use perturbation theory in several different ways. In general we partition 
the Hamiltonian in eq. (2.8.1) into 

Je = Je Q + , (2.8.46) 

where we can solve the Schrodinger equation for and we regard J€\ as a perturbation. 
The first application area is presented in section 2.9 when is the molecular Hamiltonian 


Download free eBooks at bookboon.com 


29 










ATOMISTIC MODELS 


MOLECULAR QUANTUM MECHANICS 


in eq. (2.4.2) and is the interaction with an electric field, i.e. the perturbation in 
eq. (2.8.46) is a physical perturbation of the system. 

Secondly, if we are not able to solve the Schrodinger equation for the Hamiltonian but 
instead we have a solution for an approximate Hamiltonian ,i£o> we can regard the difference 
between the exact and the approximate Hamiltonian as a perturbation, 

Je x =Je-JeQ, ( 2 . 8 . 47 ) 

and apply perturbation theory to include corrections to J&q. The typical example is that we 
cannot solve the Schrodinger equation for the molecular Hamiltonian in eq. (2.4.2). Rather 
we introduce the Hartree-Fock approximation as in section 2 . 10 , and in section 2.12.2 we 
use perturbation theory to correct for the Hartree-Fock approximation in Moller-Plesset 
perturbation theory. 

Finally in a very similar use case, we do a Taylor expansion of the Hamiltonian around JCq, 

00 y 

= — x n Je n , (2.8.48) 

n^l n - 

where we again assume that we can solve the Schrodinger equation for and we include 
J6 n as perturbations order by order. Here the standard example is vibrational motion 
in molecules where we do a Taylor example of the potential energy surface, and is 
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the harmonic oscillator, given for a diatomic molecule in appendix 2.A.2, and the leading 
correction from perturbation theory arises from the anharmonicity of the potential energy 
surface. 


2.9 First and second-order electric properties 


2.9.1 Multipole expansion 


The starting point is classical electrostatics where we regard a set of point charges interacting 
with a test charge. The interaction energy, V, between a charge distribution described by a 
set of N point charges, cfr , and a test charge, q r is given by Coulomb’s law as (using atomic 
units) 


N 

v=L 

i =1 


didt 
-R— 


(2.9.1) 


it 


where Ru is the distance between particles i and t. Note that the Coulomb interactions 
within the set of N point charges are not included in eq. (2.9.1). The electrostatic potential 
in a point t, <p t , is defined from the interaction energy V obtained by placing a test charge q t 
in point t, 


V = q t <Pt > 


(2.9.2) 


which leads to the following definition of the electrostatic potential in point t from the 
charge distribution, 

N a 

<pt = Zir ( 2 . 9 . 3 ) 

i =l R it 

Consequently, the electrostatic potential at point charge i, (pt, arising from the test charge, 
t, is given as 

(pi = (2.9.4) 

Kit 

so that the potential energy can be written as 

N 

V=Y J (pidi■ (2-9.5) 

i=1 


The next step is to carry out a multipole expansion, i.e. a Taylor expansion of the electrostatic 
potential at each point i around a common origin where it is assumed that | Roi I « \Rto\ for 
the Taylor expansion to converge rapidly. For a Taylor expansion in Cartesian coordinates, 
the electrostatic potential at particle i,tpi, becomes (see figure 2.5) 


_ dt _ dt 

\Ru\ I Rto + Roi\ 

l 

— (pO + Ri,afil i,a(PO + ^ Ri,<xRi,(fl /,/i ^ i,a(PO + • • • 
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Figure 2.5: Molecular multipole expansion. For each particle (electron or nucleus) i, a Taylor expansion of the 
electrostatic potential is carried out around a common origin Rq . Note that Ri is used to denote a point so that 
the distance vector is given by (^/ - Ro^- In general, the explicit reference to the origin is however dropped 

and thus Ri may be interpreted both as the point and the distance vector. The explicit reference to the origin 
is only needed in a few cases, for example to investigate the origin dependence of molecular properties (see 
exercises 2.7 and 2.8). 


where cpo (also written as (p ^ } ) is the electrostatic potential evaluated at the origin Rq given 
as q t l\R t o and a is the mth derivative of the electrostatic potential with respect 

to a Cartesian coordinate x, y or z calculated at the origin O. The Einstein summation 
convention is used for the subscripts a and /3 denoting the Cartesian coordinates x, y and 
z. In the last line of eq. (2.9.6), a special notation is used for m = 0, where the whole term 
is reduced to cp^. Using the Taylor expansion of the electrostatic potential in eq. (2.9.6), the 
interaction energy, V, in eq. (2.9.1) becomes 


N 


V ~ ^ tfityO + tfiRiyttVO'Cc p + * * * 


i=1 
N 


f N \ N \ i N \ 

L^ko+ JL^i R ha\(Po] a + - JLyi R Ua R i,p\<Po ) a p + 
\i= 1 ) Vi=l ) Z U=1 ) 


_ mol mui uj , - rjiiiui , 

- q (p 0 + q a <Po,a + 2 ^O.a/3 + ' 

oo 

= L 


,mol,„(l) 


^mol,_(2) 


m =0 


ml 


M [m) wW 

1Y1 ai...a m To,ai...a n 


(2.9.7) 


where an electric moment of order m, is defined as 


M« s 

UL \ • • • LI ft 


N N m 

= i =e n«<.« 

;=1 i= 1 p=l 


(2.9.8) 


In particular, the charge, q (m = 0), the dipole moment, q a (m = 1), and the second moment, 
Q aj g (m = 2), are defined as 


N N N 

q = y (fi > E-a = y qiRi,a > Qa(i = y • (2.9.9) 

i= 1 i=l i=l 

The electric field, E^ a , in point k is defined as minus the gradient of the electrostatic 
potential in the same point k, 

E k,a = -<p [ ^ a , (2.9.10) 
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which is consistent with the definition of the force as F a = - V a V, so that the electric field is 
the force acting on a particle with unit charge. The electric field gradient, Ek, a p> is trivially 
defined as 

E t . al> = -«4V (2.9.11) 

and in general the m-th derivative of the electric field, Ek >ai ...a m > is defined as 


Ek,a\...a m < Ek,a\...a m 

The interaction energy in eq. (2.9.7) is written in terms of the electric field as 


(2.9.12) 


1 


V WO ^ai...a„ 


dm) 


jpim) 

0,oc\...cc m 


(2.9.13) 


The electrostatic potential at the origin, (po, from a point charge, q t , is given by using 
eq. (2.9.4) as 

<,Po = ^ L = qtT {0 \R Ot ), (2.9.14) 

Eot 

where the T -tensor (or interaction tensor) of rank 0 is defined as 

T {0 \Rij) = 2_. (2.9.15) 

R ij 
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We now introduce the more compact notation, 

T (0) (Rij) = Tf) . (2.9.16) 

The gradients of the potential as well as the electric field and its gradients are thus obtained 
from the gradients of the T -tensor. The T -tensor of rank 1 is defined as 


T a) 

1 ij,a 


= y • T l 

- V 7,« 1 ;: 


(0) 


i] 


^j.a t> 


R 


l J 


Rij.a 



(2.9.17) 


where Rij, a is a component of the vector connecting particles i and j, Rij >a = Rj, a ~ Ri,a- 


Example 2.1: Force on a point charge 

Since the potential energy V for a charge cj j in an electrostatic potential (pj 
eq. (2.9.5) as V = qj(Pj, the force Fj iU is given as 

is given by 

II 

$ 

> 

1 

II 

S 

tC 

(2.9.18) 

where we in the second step have used the definition of the electric field in eq. (2.9.10), 
the electric field is the force on a point charge with the value +1. If the electrostatic 
potential arises from a charge q. 

— S 

E-r 

1 

II 

O 

hT 

s 

> 

1 

II 

(2.9.19) 

so that 

Ej, a = -QiTifa 

(2.9.20) 


Since R ijt a = -Rji.a, we have 


T m _ _t’(d 

ji, a 1 ij,a' 


(2.9.21) 


The T -tensor of rank 2 is consequently given as 


r r{ 2) 

ij.uP 


= V 


j.p 


J’(l) _ 

ij.a ~ 


3Rij,aRij,p fiapR'ij 



(2.9.22) 


where 8 aj g is the Kroenecker delta function. The T (2) -tensor is traceless, 


j^(2) _ J,{2) _|_ j(2) _|_ ji( 2) _ q 

1 ij f aa ~ 1 ij,xx ' 1 ij,yy 1 ij,zz ~ 


We also have 


2^(3) 

ij>aPr 



1 5Rij } aRij } l3Ri j 




Ri j,afi fiy + Rijyfifiay + Ri j y y^ap 


(2.9.23) 


(2.9.24) 
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and 


t (4) 

ij.a/iyS 


r.9 fl Q5Rij,aRij,pRij,yRij,S 


R iJ K 


1 |^zj,a^z7,/^7<5 + Rij,aRij,yfip5 + Rij } aRij } 8$ft 7 

+ ^z j,pRijyj^a8 + ^z j,pRij } 8^aj + ^z j t y^ij,8^ap j 

+ 37?l(<VV + Sayfi05 + ^aS^/ly) j • 


In general, we have 




ij,a\...a n - V ;',a»--- V y,ai 


Eqs. ( 2 . 9 . 21 ) and eq. (2.9.23) are generalized to 


ji{n) 


= (-1 r r" 

ji,a\...a n v ' z;,ai...a„ 


" ' r( " ) and 


1%. 


rj-'ill) 

ij,ai...a m a m ...a n 


= 0 


(2.9.25) 


(2.9.26) 


(2.9.27) 


Normally, the second moment in eq. (2.9.9) is instead given as its traceless counterpart, the 
quadrupole moment. As an illustration, the interaction energy, V, between a test charge, c/ r , 
and the second moment of the charge distribution is regarded, 


V - 2 Q°> a P■ 

A constant contribution, A, is added to each of the diagonal terms of Qo, a p> 

Qo,aP = Qo,af) + , 


(2.9.28) 


(2.9.29) 


The modified interaction energy, V, becomes 

V = Tolaf,^ = V + i A = V • 

where we have used that the T (2) -tensor is traceless (eq. (2.9.23)). It is thus noted that the 
trace of the second moment does not contribute to the interaction energy since A can be 
chosen as 

A = “ Q rr . (2.9.31) 

Therefore, the quadrupole moment is defined as (Buckingham 1967) 

®a/3 = 2 Qafi ~ 2 aft > (2.9.32) 

where a factor of 3/2 is introduced to be consistent with the literature. Since the quadrupole 
moment is traceless, it only contains five independent components. This is equivalent 
to d-orbitals for describing atomic orbitals, where we have five orbitals in spherical 
polar coordinates as compared to six d-orbitals in Cartesian coordinates. Similarly, the 
quadrupole moment describes a second-order shape/anisotropy of the charge distribution, 
and the term x 2 + y 2 + z 2 is isotropic. 
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2.9.1.1 Units 

Even if SI units are recommended, they are often not suitable in chemistry since they become 
inconviently large or small. A typical length of a covalent bond is in the order of 0.1 nm, 
which indeed would be a suitable unit, but historically we instead use angstrom (1 A = 
0.1 nm) so a typical bond length is around 1-2 A. In quantum chemistry, we also commonly 
use the atomic unit system, which was introduced in table 2.1, where for example the charge 
of an electron is -1. 

For dipole moments, we normally use the unit debye, which is a CGS unit and also has 
its origin in the electrostatic unit (esu) system. 1 D is the dipole moment of two charges 
with different sign, but equal magnitude of 0.2081943 e, separated by 1 A. The conversion 
factor from atomic units thus becomes 1 e • bohr = 2.541766 D. In SI units, we have 1 D = 
3.33564 • 10 -30 C-m. A typical molecular dipole moment is of the magnitude 1 -10 D. Using 
the smallest available SI prefix yocto for 1 • 10 -24 , would still mean that the magnitude of a 
typical molecular dipole moment would be around 1 • 10 -6 yC-m, illustrating what we above 
meant with “inconviently small”. 

For quadrupole moments, the buckingham unit is used, where 1B = 1 DA, i.e. the quadrupole 
moment of two dipoles of the magnitude 1D, but opposing directions, placed with a distance 
1 A from each other. The conversion factor from atomic units becomes 1 e • bohr 2 = 
1.344911 B. 
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Ro 


Figure 2.6: Molecular multipole expansions for intermolecular interactions. The distance vector R ms is defined 
as Rah = Rq,b ~ Rq,a■ Furthermore, Rj = fj - Rq and Rj = r,- - Rq. 


2.9.1.2 Multipole expansion of two interacting molecules 


Two charge distributions, A and B, are considered, where system A has Na point charges and 
system B has Nb point charges. The Coulomb interaction energy between the two charge 
distributions is 


n a n b 

Vab = E E 




=ij=i \Rj-Ri\ 


i=l j=\ \rj - Ti+R A B\ 


(2.9.33) 


where the distance vectors are defined in figure 2.6. If fj and fj are small compared to Rab> a 
multipole expansion may be carried out around fj = 0 and fj = 0. The multipole expansion 
of two interacting charge distributions becomes 


Vab 


g g R_iRj 

hi h \ fj-n + R\ 

Na Nb r j 

= L Zk' T ABRj + Ri T AB,aRj r j’ a + nR‘ T AB,apRj r i> ar j’P + ''' 

- qm,a TaIccRJ ~ Ri r i,a ^ab^R A j,(i ~ ^RAi,a T^^Rjfj.pfj.y + ■■■ 
+ ~ Ri r i,a t ir p T^ a p qj + \qiri, a ri,pT™ iafyr qjrj, r 
+ \Riri,ari,pTf Baprd qjrj, r rj, 6 + ... J , 


(2.9.34) 


where the definition of the T -tensor in eq. (2.9.21) has been adopted. It is noted that minus 
sign appearing in the linear terms in rj ta (and subsequently also for other odd order terms 
in r J)a ) arises from the definition of Rab = Ro,b ~ Ro,a as pointing from Ro,a to Ro,b (see 
figure 2.6). Defining an electric moment of order m as in eq. (2.9.8), eq. (2.9.34) can be 
rewritten as, 


N a N b 

Vab = E E — 
i=U=i I r i 


RiRj 

-n + R\ 


oo 

= E 

m,n =0 


——— 

i } 0C\...0Cyyi AB,cc\... ttm+n j yttm+l-’-ttm+n 


(2.9.35) 
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2.9.2 Polarizability in tin external field 

The dipole-dipole polarizability, a a p, is defined as the linear response to an electric field, 
Ep, (Buckingham 1967) 

Ha d = a ccpEp , (2.9.36) 

where /j" ld is the induced dipole moment. Frequently, we use a model of representing a 
molecular systems by a set of point particles, and we denote the polarizability for particle i 
as oci tCC p. Analogously to eq. (2.9.36), we define cti }(X p by 

p‘ nd = a ita pE i>p , (2.9.37) 

where /il nd is the induced dipole moment of particle i and E-^p is the local electric field 
at particle i. The electrostatic interaction energy, V e \ e , between a dipole moment and an 
electric field is given in eq. (2.9.7) as 

Tele = ~/jp d Ep • (2.9.38) 

In classical theory, there is, however, an additional energy term arising from the work 
required to create the induced dipole moment, the self-energy V^if. The polarizability 
describes the mobility of charges within a charge distribution. The work V^if required to 
move a charge, q, in an external force, Fp, is 

R 

V sel f = f FpdRp. (2.9.39) 

o 

Using that the force is Fp = qEp and the induced dipole moments is dq'^ d = qdRp gives 

T^.ind rii n d 

r* r 1 

Vseii = / Epd^ B ' 9 = 36 > j („„„= = (2.9.40) 

0 0 

where we have used partial integration in the second last step (see exercise 2.9). The total 
energy termed the induction energy, V^, is the sum of the electrostatic energy and the self¬ 
energy, 

kind = -R l p d Ep + \l^p d Ep = -\^p d Ep = - l -a a pEpE a . (2.9.41) 

The self-energy thus cancels exactly half of the electrostatic energy. It is furthermore noted 
that the induction energy is quadratic in the electric field. 

2.9.3 A molecule in an external potential 

2.9.3.1 Dipole and quadrupole moment 

To study a molecule in an external electrostatic potential, Rayleigh-Schrodinger perturba¬ 
tion theory (see section 2.8.1) is adopted. The Hamiltonian for a molecule in an external 
potential may be written as 

Je = Je mQ \ + A f , (2.9.42) 
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where jfe mo i is the molecular Hamiltonian discussed in section 2.2, A is an order parameter 
that determines the order in the perturbation expansion, and Y is the interaction operator 
between a molecular charge distribution and an external electrostatic potential, 


N 


v = E Zjt P ] ~YjYi > 

1=1 i=1 


(2.9.43) 


where the minus sign arises from the charge of the electron. As in section 2.2, we regard 
a molecule consisting of N nuclei and n electrons. Using the multipole expansion of the 
electrostatic potential in eq. (2.9.6), the interaction operator in eq. (2.9.43) becomes 


/ iV n 

^ = E Z lRl,ai ■ ■ ■ Rl,a m ~ E n,ai • ■ ■ r i,a m ] Yq 


m =0 \/=l 


im) 

ai...a n 


(2.9.44) 


The operator for a molecular electronic moment is identified as 

N m n m 

^S..a m = E Zl n R I,ap - E n n,a p • (2.9.45) 

1=1 p=l i=lp=l 

Consequently, the first-order energy in eq. (2.8.9) using Rayleigh-Schrodinger perturbation 
theory becomes 

CO 

(2.9.46) 


4 u =<o*vr> = 1 

m=0 m - 


By comparison to eq. (2.9.7), the molecular charge, q mo \ (m = 0) is trivially given as 


„mol 


N 


= Y J Z I -n, 
1=1 


the molecular dipole moment, /i™ 1 , (m = 1) is identified as 


N n 

.mol _ r 7 d , a.„(0 ). 


Ha = E Z l R ha + (To I - E n,a Wo ) > 

1=1 i=1 

and the molecular second moment, Q™ 1 , (m = 2) is identified as 


(2.9.47) 


(2.9.48) 


N 


Qap = E ZiRi,aRi t p + << } l - £ r Ua r i} p |^ 0) ), 


1=1 


i=1 


The molecular moment of order m is given as 


N 


M nU„ = E z l R lm ■ ■ •+ (V'S 0 1 - E -r „ m !<>). 


(2.9.49) 


(2.9.50) 


i=i 


i=1 


The first-order correction to the energy becomes 


_ mol , mol (1) , x n mol (2) , 

- q (po + Ha ( P 0t a + 2 Q ^ ( P 0ta p + --- 

_ ^mol tn , ,mol c ^ /^mol r? , 

— q cpo [i a Eo,a + 


co 1 

= E 


M [m) (D [ ™ ] 

iV±rv ' "-rO,ai...a„ 


- yyj] oc\...a m y r 0,a\ 
m =0 rn - 


(2.9.51) 


It has thus been demonstrated by comparing to eq. (2.9.7) that the first-order energy 
corresponds to classical electrostatics. 


Download free eBooks at bookboon.com 


39 



ATOMISTIC MODELS 


MOLECULAR QUANTUM MECHANICS 


2.9.3.2 Electrostatic potential 


Rewriting eq. (2.9.43) as the electrostatic interaction between a molecule and a test charge 
q t , we get 


f= yZiq 1 £ q t 
h R n h r it ' 


( 2 . 9 . 52 ) 


The electrostatic potential at the test charge, (p t , is given by eq. (2.9.3) leading to 


. £ Z/ f 1 

<Pt=Ly-~L — ’ 

i=i K it i=i r n 


( 2 . 9 . 53 ) 


so that 

f = <p t qt • ( 2 . 9 . 54 ) 

Again using the first-order correction to the energy from Rayleigh-Schrodinger perturbation 
theory in eq (2.8.9), 

4 1] = {yfWtyf) = tyT\Qt\vf)qt , ( 2 . 9 . 55 ) 

where thus the electrostatic potential in point t from a molecule becomes 


V, = <vS”i <Mv4°’> = £#--<vS»i 1 -V«»>. 

/=1 tilt i =1 r it 


( 2 . 9 . 56 ) 
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2.9.3.3 Polarizability 

For the second-order contribution to the energy in eq. (2.8.18), 


1 


p=l£ 0 -E p 


(2.9.57) 


it is noted that only the terms in T including explicit dependence on the electronic 
coordinates, r*, a , are non-zero since the electronic states are orthogonal, i.e. (y/^= 

*>pq 


lUllldlCOj / i f (X) ClIC olllLC L11C C1CL/L1 UlllU oldlCo die Lll lllUgUlldlj l.C. \l £ p \ Cf 

vpq. Thus the term for m = 0 in eq. (2.9.44) as well as all nuclear contributions to 
eq. (2.9.45) vanish. The second-order energy becomes 


in 


OO OO I 

4 2) =LL 


m 


= 1^1 mini 


CO | 

£ 15 )—m <■< <>>« i-CV. ^ 0) > 

P=l£ 0 -Ep 


X m [m) <n {n) 

VO,a l ...a m Vo,p l ...p n ■ 


(2.9.58) 


Adopting the following notation for a polarizability, P 


r {m } ri) 

ai...a m ,Pi...p n ’ 


CO 1 

= 115— s (V'S” \-d™..a m l< > XV'™ 1-4)”’. IrS”) 

p=\£ n ~£ n 


D (m,n) 

CC\...(Xm>Pl-“Pn c (0) _v)- 

p=l£ 0 -£ p 


(2.9.59) 


eq. (2.9.58)1 is rewritten as 


CO CO 1 

p (2) _ V- y' 1 p{m,n ) Am) An) 

o o f oc\...oc m o,/3i.../3 w 

m= 1 n= 1 ,n - n - 

The leading term (m = 1, n - 1) becomes 

\ 

—7iT^() 0) | - E '7,al'l/ / p ) >^o, ) a (V / p ) l - E ■ (2.9.61) 

£p i =1 i =l ’ 

Defining the molecular dipole-dipole polarizability (m = 1, n = 1), a™ 1 , as 

1=1 1=1 


(2.9.60) 


CO 

4 2) = I 

P = 1 


<C 1 = 2 £ 


P=1 


P (0) _ AO) 
c p c 0 


(2.9.62) 


where we note the factor of 2 and the reverse of sign in the denominator as compared to 
eq. (2.9.58). This leads to 


£ 


( 2 ) 

0 


_i mol (1) (1) , 

^0,^0,a + *** 

~ 2 Eo,pEo,a + • • • > 


(2.9.63) 


so that the expression for the induction energy in eq. (2.9.41) has been obtained. This 
second-order energy quadratic in the electric field is therefore termed the induction energy. 
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2.9.4 Frequency-dependent polarizabilities 

We use the general result from time-dependent perturbation theory in eq. (2.8.34), 

t 

\f\yf) J /(f')e iWmpr 'df', (2.9.64) 

o 

so that the first-order correction to the time-dependent wavefunction is 

v P (1) (f) = X4 1) ( f )e _iWfcf l^° ) ) • (2.9.65) 

k 

The Hamiltonian, is as in eq. (2.8.19), divided into 

Je{t ) = + XJ€ x {t ), (2.9.66) 

where (t) - T fit). Here we regard the response to an oscillating electric field, 

=2f e Et cos(Mt) =f (e {£+i(0]t + e {£ - i(0]t } , (2.9.67) 

where the term e Et ensures that the field is turned on in a distant path. The parameter e is 
supposed to be small and will be allowed to approach zero giving a steady-state response 
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when the field has been turned on a for a long time. Now we intgrate eq. (2.9.64), 


a^\t) 


X 

= J e {E+i0)k P +i(0)t ' + e (e+iw *p- iw)f ' df 


_1 / e {£+\(x) kp +\(x))t' e {£+\(x) kp -\(x))t'' 
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<vri y\rp 
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(yf\r\Vp) 
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-— + -— I 


£—•0 
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(Olc p + (x)-\£ (Ojcp-iO-\£ 
f 0 i (ujcp +(*>) t e i ico kp -( d) t \ 


+ 


G)jcp + (O 0)] C p (O 


We now regard time-dependent dipole moment of the ground state (p = 0), 

vaW = mmPamt)), 

with 

oo 

Vit) =y/ { Q ] e~ i(0ot + £ Aa [ ve~ il ° kt 

lc^ 0 

Considering each order i in A*, we get for i = 0, 

A*!? = 

and for i = 1, 

OO , \ ^ 

a 4 1} = £K 1)(f) 

1c? 0 

Using eq. (2.9.68) with T - -ppEp and p = 0, 
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Using the definition of the polarizability in eq. (2.9.36) slightly extended to 
frequency- dependence, 


= OC a 


p{(o)Ep | 


e icot + e -icot 


)■ 


leads to the following result for the frequency-dependent polarizability, 
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2.10 The Hartree-Fock approximation 


The starting point for the Hartree-Fock approximation (Hartree 1928, Fock 1930, Roothan 
1951) is the energy for the Slater determinant in eq. (2.6.36) using the more compact notation 
in eq. (2.6.34), 

Eq = Z hj + — Z f Jij — Kij | + V nn 
i =1 2 i,j =1 

n „ 1 n , x 

= + - Z [(XiXj\g\XiXj)-(XiXj\g\XjXi)] + Vnn- (2.10.1) 

i =1 Z i.j'=l 

The Slater determinant is an approximate wavefunction and we can thus minimize the 
energy in eq. (2.10.1) using the variational principle in section 2.7 by optimizing the 
orbitals. We would, however, like to retain that the orbitals form an orthonormal set, and 
this condition is imposed by a constrained minimization using the method of Lagrangian 
multipliers. The Lagrangian, Eq, becomes 

Eq = E 0 - Z Aij^(Xi\Xj)—dij^=EQ— Z A j j ^ S j j — 8 j j J , (2.10.2) 

i.j =i 

where we have one Lagrangian multiplier Afor each pair of orbitals and Sjj is the 
Kroenecker delta function indicating the required orthonormality of the orbitals. We define 
an element of the overlap matrix S,y as 


$ij ~ (Xi\Xj) ■ 

A small variation of the Lagrangian, 8Eq becomes 

8Eq = 8Eq- £ A ij[{8xi\Xj) + <Xi\SXj>) - 
i,i =1 

where a small variation of the energy Eq in eq. (2.10.1) becomes 


(2.10.3) 


(2.10.4) 


8E 0 


Y,(Sxi\h\Xi) + (Xi\h\8xi) 

i =1 

\ Z (( s XiXj\g\XiXj) + (Xi$Xj\g\XiXj) + (XiXj\g\8XiXj) + (XiXj\g\Xi8Xj)) 
Z id=l 

(^Aa;lglAyA/> + (XiSXj\g\XjXj) + (XiXj\g\8XjXi) + <XiXj\$\Xj8Xi>) 
Y,( s Xi\h\Xi) + (Xi\h\8xi) 


i =1 
n 


+ 


Z - feAjlglA/A*) - to/lgl^OT) (2.10.5) 

*J=1 


In the last step, we have used that i and j are just dummy indexes and can be interchanged. 
Introducing two operators, the Coulomb operator, 


Jj\Xi> = <Xj\g\Xj)\Xi> > 


( 2 . 10 . 6 ) 
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and the exchange operator, 

rfjlXi) = (Xj\g\Xi)\Xj) - (2.10.7) 

where the operator exchanges the orbital it operates on. We thus get 

6E 0 = t,(8xi\h\Xi) + (Xi\h\8 X i) 

i =1 

+ X ~ (Xi~ \S%i) • (2.10.8) 

hj = 1 

Defining the Fock operator / as 


/=^+E( < /7-^) - (2-10.9) 

,/= ] 

so that the total Hartree-Fock Hamiltonian, ^hf» becomes 

^HF =£,fi, (2.10.10) 

i =1 

where the subscript i denotes that we have one Fock operator for each electron in the system. 
We get 

8E 0 = it(6xi\fi\Xi) + < Xi\fi\8Xi > (2.10.11) 

i =1 

so that the total Lagrangian becomes 

5£ 0 = E<^l/dB> + <ld/d^>- E Xij 

i =1 ij=l 

Assuming that either a small variation or leads to that the variational principle is 

fulfilled, i.e. 8Eq = 0, gives two relations to be fulfilled simultaneously, 

L(8xi\fi\Xi) - E *ij<8Xi\Xj) = 0 - (2.10.13) 

i =1 i,j =1 

and 

f = (2.10.14) 

i =1 *d=l 

Using that 

(2.10.15) 

and subtracting eq. (2.10.13) from the complex conjugate of eq. (2.10.14) from each other, we 
get 

E (a / ;-A; i .)<^ / |^> = 0, (2.10.16) 

hj =1 

so that the condition is fulfilled if Ais an element of a Hermitian matrix. In the next step, 
eq. (2.10.13) is rewritten as a set of eigenvalue problems rather than an expectation value, 


[(8Xi\Xj) + <Xi\8Xj)) ■ (2.10.12) 
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fi\Xi) = L*ij\Xj) V/. (2.10.17) 

7=1 

We can now make a specific choice of the orbitals in eq. (2.10.17) by a unitary transformation 
so that A ij becomes a diagonal matrix. With the notation c,- = Xu, 

fi\Xi)=ei\Xi) Vi, (2.10.18) 

where this choice of orbitals is termed the canonical orbitals and ey is an orbital energy. 
Eq. (2.10.18) is referred to as the Hartree-Fock equations and is a set of coupled eigenvalue 
equations since the Fock operator in eq. (2.10.9) depends on all the orbitals through the 
Coulomb and exchange operators in eqs. (2.10.6) and (2.10.7), respectively. An iterative 
method, the self-consistent field (SCF) method is therefore commonly employed with a set 
of starting orbitals to get an initial guess of Eq. (2.10.18) is then solved for this set of /,• 
providing an updated set of orbitals and thereby also an updated and the procedure is 
repeated until convergence is reached. 

Rewriting the quantum mechanical expression for the molecular electrostatic potential in 
eq. (2.9.56) in terms of a Slater determinant, 

N Zt n 1 

<pt = Y.ir ~L<Xi\—\Xi) , (2.10.19) 

1=1 K It i =i Ht 
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where the second term on the right-hand side is the integral that appears for the Coulomb 
operator in eq. (2.10.6). The averaging over one of the electronic coordinates in the 
definitions of the Coulomb and exchange operators in eqs. ( 2 . 10 . 6 ) and ( 2 . 10 . 7 ) thus leads 
to that the two-electron Coulomb term in eq. ( 2 . 10 . 5 ) has been replaced by the interaction 
between one electron and the average electrostatic potential in eq. ( 2 . 10 . 8 ). We refer to this 
kind of approximation as a mean-field approximation. 

The orbital energy e* can be obtained as an expectation value, 

C = (Xi\fi\Xi) ~ hi ^ \jij ~ ^ , ./| ’ (2.10.20) 

7=1 

where we instead of integrating over the Fock operator /,■ return to the initial two-electron 
integrals. Consequently the energy for a Slater determinant in eq. (2.10.1) becomes within 
the Hartree-Fock approximation 


Eq — y. £j _ ^ {jij Kij\ + V nn , (2.10.21) 

i =1 Z i,j =1 

and is thus not just the sum over the orbital energies. In eq. (2.10.21), we have also added the 
classical nucleus-nucleus interaction energy. 


2.11 Basis set expansion 

We have to represent the spin-orbitals mathematically in some way to be able to do actual 
calculations to solve the Hartree-Fock equations in eq. ( 2 . 10 .18). The set of basis functions, 
in quantum chemistry termed the basis set can be chosen in different ways. We will return to 
different choices for the basis set, here we just assume that a spin-orbital, %i is expanded in 
m basis functions <p p (Roothan 1951), 


m 

I Xi) ~ Y Cjv\(pp) > (2.11.1) 

P= 1 

where Ci P is an orbital coefficient. The Hartree-Fock equations in eq. (2.10.18) become 

m m 

fi ^ Cip\(pp) —Gi ^ Cip\(pp) V i . (2.11.2) 

P =1 7=1 

We introduce the following matrix notation 




<P= (01,02. • 

• • 0m) 

y 


(2.11.3) 


r C\i 


( 

Cll 

Cl2 •• 

\ 

c l n 


c i ~ 

C2i 

and C = 

C21 

C22 

c 2n 

(2.11.4) 


y Cmi j 


K Cml 

c m2 

Cmn 
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so that we can write 




>< 

II 

o 

i and 

rfc. 

We also define the Fock matrix as 




r <0iiM> 

<0ll/l02> 

... <0ll/l0m> 

F = 

<02l/l01> 

<02l/l02> 

... (Cp2\f\cpm) 


k (<Pm\f\<Pl) 

<0ml/l02> 


In addition, the overlap matrix is defined as 




r <0ll01> 

<01102> 

• •• (0110m) 

S 


<02101 > 

<02l02> 

(0210m) 




<0ml02> 

(0ml0m) J 


(2.11.5) 


( 2 . 11 . 6 ) 


(2.11.7) 


Again applying the variational principle and the Rayleigh-Ritz method in section 2.7, we 
arrive at (see exercise 2 . 10 ) the Roothaan-Hall equations (Roothan 1951, Hall 1951) 


Fc; = CjSci or FC = SCe, 


( 2 . 11 . 8 ) 


where e is a diagonal matrix with the orbital energies e,-. We have thus converted the Hartree- 
Fock equations in eq. ( 2 . 10 .I 8 ), which are a set of coupled integro-differential equations, 
to finding the eigenvalues and eigenvectors of the Fock matrix. The overlap matrix S is 
regarded as the metric of the Fock matrix F since it reduces to the unity matrix 1 for an 
orthonormal basis. As for the Hartree-Fock equations in eq. ( 2 . 10 .18), the Roothaan-Hall 
equation in eq. (2.11.8) is solved by an iterative method since the Fock matrix F depends on 
the eigenvectors, i.e. the orbital coefficients, in C, which we also here refer to as the self- 
consistent field (SCF) method. 


2.11.1 Density matrices 

To introduce density matrices, we regard an element of the Fock matrix F pq in eq. (2.11.6), 
Fpq = <' <Pp\f\<Pq) 

71 

= {(pp\h\(pq) + y, ((pplcflj ~/ I0q) 

7=1 

n 

= {(pp\h\(p q ) +Y,( ( PpXj\g\<PqXj) ~ < (p P Xj\g\Xj(pq > 

7=1 

n m . 

— ((pplhlcpq) + y y Cj r Cj s [{(Pp(p r \g\(pq(p s ) — {(Pp(p r \g\(ps(pq) 

j=l r,s=l 

TYl 

— {(pp\h\(pq) + D rs i{(pp(f) r \g\(f)q(p s ) — {(pp(p r \g\(p s (pq)\ y (2.11.9) 

r,s= 1 
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where we have defined the density matrix as 

n 

Drs = CjrCjs • (2.11.10) 

7-1 

We also introduce the following notation for the two-electron integrals, 

Gprqs = ((pptprlsltpqtps) ~ i^Pp^rlgl^Ps^Pq') ( 2 . 11 . 11 ) 

so that 

m 

Epq = hpq + DfsGprqs • ( 2 . 11 . 12 ) 

r,s =1 

The energy for a Slater determinant in eq. ( 2 . 10 . 1 ) is recalled and rephrased in terms of 
density matrices, 


Eo 


n ^ 1 n 

Y.(Xi\h\Xi) + ~ E \(XiXj\g\XiXj) - (XiXj\g\XjXi)\ 

i =1 Z *,7=1 

n m \ n 171 

L L c ip Ciq(*Pp\fl\(pq) + - 9 L L c ip c iq c jr c js 

i-lp } q-\ ^ i,j=l p,q,r,s=l 

m y m 

Y Dpqhpq + — Yj DpqD rs Gp r q S + V nn . 
p y q =1 ^ p,q,r,s=l 


+ Vnn 

(< (pp<Pr\g\<f>q(ps) ~ < (pp(pr\g\<Ps<Pq >) 

(2.11.13) 
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2.11.2 Basis sets 

We need to have a mathematical representation of the basis functions <p p {i) in eq. ( 2 . 11 . 1 ). 
Traditionally there are two approaches, either plane waves are used which is suitable when 
applying periodic boundary conditions in solid-state chemistry or as a linear combination 
of atomic orbitals (although these “atomic orbitals” may not be a solution to the Schrodinger 
equation in eq. ( 2 . 5 . 2 ).). Two types of atomic orbitals are commonly employed, Slater-type 
orbitals (STOs) that decay exponentially with the distance from the nuclei, and Gaussian- 
type orbitals (GTOs) which adopt Gaussian functions. 

The choice of GTOs versus STOs as the functions in the basis set is still the subject of a vivic 
discussion, although GTOs dominate today in the commonly employed quantum chemical 
software. An STO is indeed the exact solution to the one-electron problem as the hydrogen 
atom, but on the other hand, we do not know how the exact solution looks like for a many- 
electron problem. STOs have a more long-range decay than GTOs, and describes the cusp 
at the nucleus. GTOs were originally introduced to describe one STO by several GTOs (Boys 
1950) since the two-electron integrals in eq. ( 2 . 11 . 11 ) are much easier to calculate by using 
GTOs as compared to STOs. 

A useful relation for a Gaussian function positioned at nucleus I, <pi,i(r), 

(pi,dr) = (^) 4 e - a '( f -*/) 2 (2.11.14) 

is its product rule (see exercise 2.11 for a derivation), 

(pi, I (p jJ = Ce-^ K]Z (2.11.15) 


where 


P = ai + a h , 


- aiRi + ctjRr 

R k = ——-— and 

ci[ + ai 


C = 



e «*■+«/ 


t Ri-Rj ) 2 


(2.11.16) 


which simplifies the integral calculations greatly for GTOs. 

In general we write a Gaussian basis function in Cartesian coordinates as 


(pi{x,y,z) = x a y b z c e a ‘ r , a,b,c> 0, 


(2.11.17) 


where a, b and c are integers and an s-function thus corresponds to a+b+c = 0, a p-function 
to a + b + c=l, a d-function to a + b + c = 2, etc. 

There is a huge amount of basis sets available in the literature and accessible by simple key¬ 
words in quantum chemical software. Here we will only discuss some general principles and 
we will use the correlation-consistent basis sets (Dunning Jr. 1989, Woon and Dunning Jr. 
1993) as an example. If we first assign the number of basis functions to each atom so that 
the GTOs can contain all the electrons of the neutral atom, e.g. for C with 6 electrons we 
need two 5-functions and one p-function (which indeed is three basis functions, p x , p y and 
p z ) denoted [2slp]. We denote the GTO describing the Is electrons of C for a core function 
and the GTOs describing the 2s and 2p electrons for valence functions. A double-^ basis set 
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H-He 

Li-Ne 

Na-Ar 

cc-pVDZ 

cc-pVTZ 

cc-pVQZ 

[2slp] 

[3s2pld] 

[4s3p2dlf] 

[3s2pld] 

[4s3p2dlf] 

[5s4p3d2flg] 

[4s3p2d] 

[5s4p3dlf] 

[6s5p4d2flg] 


Table 2.3: The number of basis functions in some of the correlation-consistent basis sets, cc-pVXZ, X={D, T, Q}. 


for the C atom contains one additional s- and ^-function and in addition one d-function, 
denoted [3s2pld]. A table for the correlation-consistent basis sets, cc-pVZX, X={D, T, Q} 
including also the number of basis functions for triple-^ and quadruple-^ basis sets is given 
in table 2.3. 

In this type of hierarchy of basis sets, double-^, triple-^, etc., we add what is called 
polarization functions for each level in the hierarchy with the purpose to improve the 
description of covalent bonds. The calculated accuracy of many molecular properties 
depends strongly on the description of the covalent bonds in the molecule, where the 
molecular geometry is a typical example. It is also beneficial to have a family of basis sets 
where the hierarchy converges towards the basis set limit, so that we can control the error 
introduced in a truncated basis set expansion. The cc-pVXZ basis sets is an example of such 
a family of basis sets that converges smoothly to the basis set limit when the basis set is 
increased in a systematic way. 

Also diffuse functions are often added to a basis set to improve the description of the 
electron charge distribution far away from the nuclei, which is important for example for 
the calculation of molecular polarizabilities in eq. (2.9.62). Correlation-consistent basis 
sets with diffuse functions are denoted as “augmented” basis sets, aug-cc-pVXZ, where for 
example for C the cc-pVDZ basis set consisting of [3s2pld] are extended to [4s3p2d] for aug- 
cc-pVDZ. Further augmentation with diffuse functions are denoted doubly augmented, d- 
aug-cc-pVXZ, triply augmented, t-aug-cc-pVZX, etc. In comparison to adding polarization 
functions, no higher angular-momentum functions are added to the diffuse basis set and 
obviously the exponent a -, in eq. (2.11.17) are very different when adding polarization or 
diffuse functions. There are also other extensions of basis sets, for example in the cc- 
pCVXZ basis sets extra core functions are added to improve the description of the electron 
distribution close to the nuclei which for example is important in the calculation of the 
Fermi-contact term in spin-spin coupling constants. 


Example 2.2: Basis set limit of the HF molecule at the Hartree-Fock level. 


Hartree-Fock calculations on the HF molecule are presented in table 2.4 using the 
NWChem software (Valiev et al. 2010), where molecular properties as the dipole 
moment,quadrupole moment, polarizability as well as the equilibrium bond length 
are presented for three hierarchies of basis sets: cc-pVXZ, aug-cc-pVXZ and d-aug-cc- 
pVXZ with X={D, T, Q, 5} (Dunning Jr. 1989, Kendall et al. 1992, Woon and Dunning Jr. 
1994). Since the basis set is increased systematically, we should be able to deduce a 
suitable basis set for each property at the Hartree-Fock level. 

We define the z-axis along the dipole moment axis of the HF molecule so that 
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basis set 

E/hartree 

Hi D 

0 ZZ /B 

«xx/A 3 

a zz lA 6 

R e / A 

cc-pVDZ 

-100.019419 

1.94920 

2.03526 

0.22905 

0.59513 

0.90148 

cc-pVTZ 

-100.058021 

1.94077 

2.09725 

0.38085 

0.73062 

0.89794 

cc-pVQZ 

-100.067695 

1.93330 

2.11739 

0.48097 

0.78760 

0.89685 

cc-pV5Z 

-100.070440 

1.93189 

2.13462 

0.54619 

0.81708 

0.89685 

aug-cc-pVDZ 

-100.033474 

1.93076 

2.15627 

0.56147 

0.82766 

0.90019 

aug-cc-pVTZ 

-100.061078 

1.92502 

2.15742 

0.63241 

0.84625 

0.89912 

aug-cc-pVQZ 

-100.068568 

1.92206 

2.15772 

0.65580 

0.85106 

0.89772 

aug-cc-pV5Z 

-100.070583 

1.92163 

2.15486 

0.66118 

0.85191 

0.89595 

d-aug-cc-pVDZ 

-100.033649 

1.92112 

2.18807 

0.66093 

0.84559 

0.89992 

d-aug-cc-pVTZ 

-100.061177 

1.92134 

2.16222 

0.66619 

0.85126 

0.89916 

d-aug-cc-pVQZ 

-100.068605 

1.92128 

2.15656 

0.66575 

0.85235 

0.89727 

d-aug-cc-pV5Z 

-100.070587 

1.92153 

2.15440 

0.66512 

0.85237 

0.89694 


Table 2.4: Basis set convergence for the bond length R e , dipole moment p z , quadrupole moment Q zz , and the 
polarizability, a xx and a zz , for the HF molecule at the Hartree-Fock level. See example 2.2 for a discussion. 


ocyy = axx and Q xx = ®yy = ~\®zz (since the quadrupole moment is traceless, see 
eq. (2.9.32)). All off-diagonal tensor elements, a a p and ® a p for a ^ P, are 0 because 
of symmetry reasons. All the properties (obviously apart from R e ) are calculated at an 
experimental geometry, R e = 0.9168 A. Thus we look at the basis set dependence of a 
property separately from the basis set dependence of the molecular geometry. (It is, 
however, common to do the opposite and calculate the properties at the optimized 
geometry for each basis set, but then we mix the basis set dependence of the property 
and the molecular geometry in the same calculations.) 

The total energy, E, cannot be measured experimentally, and its inclusion in the table 
only serves as a consistency check. Since the variational principle in section 2.7 is used 
in the Hartree-Fock approximation, the energy must be lowered if the size of the basis 
set is increased. That is true both for the sequence D^T^Q^5 as well as for cc—aug- 
cc^d-aug-cc, and as can be seen in table 2.4 these conditions are fulfilled. If not, it 
would be a clear sign of that something is wrong with the calculation, and since human 
errors are prone also among computational chemists, it is a good advice to inlude 
as many consistency checks as possible. We also see for example that improving the 
basis set from aug-cc-pVTZ to d-aug-cc-pVTZ lowers the energy by only 99 /ihartree, 
whereas it has a relatively large effect on the polarizability. In general, diffuse functions 
do not decrease the total energy, E, by a lot if the basis set is well converged with 
respect to polarization functions, whereas diffuse functions may have a crucial effect 
on for example the polarizability where the tail of the orbitals needs to be described 
accurately. 

The bond length, R e , is well described by the cc-pVTZ basis set and additional diffuse 
functions have very little effect. In contrast, the polarizability is well described by the 
d-aug-cc-pVDZ basis set, and addition of extra polarization functions have a minor 
effect. The dipole and quadrupole moments are somewhere in between and in this 
example (the HF molecule using the Hartree-Fock approximation) an aug-cc-pVTZ 
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basis set gives a good description. 

Normally, we do not present results for the Hartree-Fock approximation and basis sets 
are designed to include also what we are missing in the Hartree-Fock approximation, 
namely electron correlation (see section 2.12). The requirements on the basis set will 
change depending on the computational method, the molecule and the property, and 
we will discuss more examples when we have established more accurate computa¬ 
tional methods including also electron correlation. That is also the reason why we do 
not include a comparison to experimental results in this example. We demonstrate, 
however, that it is important to have hierarchies of basis sets approaching the basis set 
limit. The actual basis set used in a study will be a compromise between accuracy and 
computational efficiency but the choice of basis set should be based on a systematic 
study of the basis set dependence. 


A large amount of various basis sets are normally included in a database within the quantum 
chemical software, and the user normally only has to provide a simple keyword to use a 
basis set in a calculation. However, many software adopting Gaussian basis sets rely on a 
common basis set exchange (BSE) server where most basis sets can be obtained (Feller 1996, 
Schuchardt et al. 2007). 

2.12 Electron correlation 

The electron correlation energy of state i, £“ rr , is pragmatically defined as (Lowdin 1955) 


jjcon _ ^exact _ ^-HF 


Tcorr 


( 2 . 12 . 1 ) 


i.e. the energy not included in the Hartree-Fock approximation, £™. The exact energy fJ® xact 
here refers to the “exact” solution of the Schrodinger equation within the molecular orbital 
ansatz for a given basis set, so the “true” solution to the problem would be to include all 
electron correlation as well as an infinitely large basis set (see figure 2 . 7 ). Since the molecular 
Hamiltonian in eq. (2.4.2) commutes with the Fock operator in eq. (2.10.9) (see exercise 2.12), 
the “exact” wavefunction t//® xact can be expressed as a linear combination of the states to the 
solution to the Hartree-Fock approximation, i/a™, 



( 2 . 12 . 2 ) 


where C JjU are expansion coefficients to be determined. Within molecular orbital theory, 
we thus focus on which Hartree-Fock states i/a™ to include and how to determine their 
respective coefficients C JjU . All methods discussed in this section are based on the restricted 
Hartree-Fock method, i.e. a closed-shell system with two electrons in each occupied orbital 
(see eq. (2.6.39)). We can therefore depict the methods by using figure 2.8 where we in the 
Hartree-Fock ground state have put two electrons in each orbital filling the orbitals with the 
lowest orbital energies (see figure 2 . 8 a). The Hartree-Fock states included in eq. ( 2 . 12 . 3 ) can 
be depicted as single excitations (see figure 2 . 8 b), double excitations (see figure 2 . 8 c), etc. 
out of the Hartree-Fock ground state. 
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Figure 2.7: Illustration of that we need both electron correlation, ideally in the full Cl limit, and very large basis 
sets to get the “true result”. 


2.12.1 Configuration interaction (Cl) methods 


In configuration interaction (Cl) methods, we rewrite the expansion in eq. (2.12.3) as 


,ci 


Vo* = CoVo +E c ^vi +••• = Qvr+EE c ?v?+E E c ?j Vfj+••• - ( 2 - 12 - 3 ) 

Id Id i a i } a, 


J>i 


b>a 


where i/a™ is the Hartree-Fock ground state (see figure 2.8a), denotes a single excitation 
of the Hartree-Fock wavefunction (see figure 2.8b) and y/ff is a double excitation of the 
Hartree-Fock wavefunction (see figure 2.8c). In the last part of eq. (2.12.3), i and j refers 
to occupied orbitals and a and b to virtual orbitals, respectively, so that yrf and yrff are 
singly and doubly excited states with the Hartree-Fock ground state as the reference state. 
Normally we truncate after single and double excitations and this method is termed CISD 
(Cl with singles and doubles excitations). The coefficients C A , are determined variationally. 
The orbital coefficients c ; y in the LCAO approximation in eq. (2.6.1) are not reoptimized, and 
therefore a truncated Cl method relies on that the Hartree-Fock approximation is a good 
starting point. 


2.12.1.1 Brillouin’s theorem 


We look at one of the energy contributions, the coupling between the Hartree-Fock state and 
all single excitations. Adopting the notation in eq. (2.10.1), 

( \ 


i a 


EE 

i a 


(Xi I h\Xa) + \ E (<XiXj\g\XaXj > - (XiXj\g\XjXa >) 


V 


7 


LL<xi\f\xa) = EE^<t*it«> = EE e ^-a = 0 - (2.12.4) 

i a i a i a 
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(c) Examples of double excitations 


Figure 2.8: Illustration of (a) the Hartree-Fock ground state, (b) single excitations, and (c) double excitations 
for a system of 4 electrons in 4 orbitals. 


where we first have used the definition of the Fock operator in eq. (2.10.9) and the Coulomb 
and exchange operators in eqs. (2.10.6) and (2.10.7), and then we have used that the 
canonical orbitals are eigenfunctions of the Fock operator / (see eq.(2.io.i8)) as well as that 
the orbitals form an orthonormal set. The relation in eq. (2.12.4) is the Brillouin’s theorem. 
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cc-pVDZ 

cc-pVTZ 

cc-pVQZ 


n 

m 

Afcci 

m 

Nfci 

m 

-ZVfci 

h 2 o 

10 

24 

6.5-10 9 

58 

8.2 • 10 19 

115 

9.4- 10 ib 

co 2 

22 

42 

9.4-10 19 

90 

9.6-10 27 

165 

1.1 • 10 34 

c 6 H 6 

42 

114 

1.4 -10 46 

264 

3.0-10 62 

510 

6.9-10 74 


Table 2.5: Examples of the number of Hartree-Fock states included in a full Cl calculation for some typical 
molecules and basis sets, n is the number of electrons in the molecule, m is the number of basis functions 
for the given basis set and molecule, and A/pci is the number of Hartree-Fock states included according to 
eq. (2.12.5). 

2.12.1.2 Full Cl 

In full Cl all possible excitations in eq. (2.12.3) are included. The number of ways, jVfci. n 
electrons can be distributed in m orbitals (where we can put two electrons in each orbital) is 
given by the binomial coefficient, 


-ZVfci = 


(2 m)! 

n\{2m- ri)\ 


(2.12.5) 


which is a factorial scaling with the size of the problem and makes full Cl a very expensive 
method also for relatively small systems (see table 2.5 for some examples). So full Cl is 
limited by both the number of electrons in the molecule as well as the size of the basis set (m 
basis functions give m orbitals). For a given basis set, however, full Cl gives the exact solution 
and is therefore valuable for validating other methods for including electron correlation. 


2.12.2 Moller-Plesset perturbation theory 

We introduced perturbation theory as a general tool in section 2.8 and here we will use it for 
including electron correlation according to Moller-Plesset (MP) perturbation theory (Moller 
and Plesset 1934). In this method, we assume that in eq. ( 2 . 8 . 1 ) is simply operator given 
for the Hartree-Fock approximation in eq. (2.10.10), 


*^HF = Lfi = L 

i=1 i=1 


hi + 


n , , 

■ E [s, 

7=1 


■x,) 


( 2 . 12 . 6 ) 


J 


where it is noted that the Coulomb, J j, and exchange, JF), terms are counted twice for each 
electron pair i and j, which is also reflected in the Hartree-Fock energy in eq. ( 2 . 10 . 21 ). The 
perturbation operator in eq. ( 2 . 8 . 1 ), J€\, becomes 


;t/ 9 _ -z/?el ' t/d _7/ 

rD7L\ — ye — tTCHp — y et 


E (Jj 

1 , 7=1 


■JfTi 


■') ’ 


(2.12.7) 


where ,i^ el is the molecular Hamiltonian in eq. ( 2 . 4 . 2 ). A rather peculiar fact about is 
that the second term on the right-hand side of eq. (2.12.7) is around twice the size of the first 
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term, whereas one in general in perturbation theory should expect the perturbation to 
be as small as possible. The zeroth-order energy, e®, becomes 

4 0) = X> (2.12.8) 

i=1 

by simply using eq. (2.10.18). The first-order correction to the energy is according to 
eq. (2.8.9) the expectation value of the perturbation operator which leads to 

= (2.12.9) 

Z *.7=1 

where the integration over the first term on the right-hand side of eq. (2.12.7) cancels exactly 
half of the integration over the second term. The Hartree-Fock energy E 0 in eq. (2.10.21) is 
thus retained as 

E 0 = £ ( ° ] +£ { Q 1] + V nn , (2.12.10) 

where we also added the classical nuclear repulsion energy. Electron correlation, as defined 
in eq. (2.12.1), thus enters first in the second-order contribution to the energy in Moller- 
Plesset perturbation theory since the sum of the zeroth-order contribution and first-order 
correction gives the Hartree-Fock energy. 

To get the second-order energy, the MP2 energy, we recapitulate the second-order energy 
from RSPT in eq. (2.8.18), 


„(2) _ V 

£ ° ~ E 7(0) 7(0) 

7=1 C 0 ~ £ i 


( 2 . 12 . 11 ) 


where we here regard a correction to the ground-state energy. As for the Cl method in 
section 2.12.1, we write the excited states, as a sum of single excitations y/f, double 

excitations etc. Because of the Slater-Condon rules discussed in section 2.6.2, higher 
excitations than double excitations do not contribute to eq. (2.12.11) for an orthogonal 
set of states. Furthermore, for the single excitations we can utilize Brillouin’s theorem in 
section 2.12.1.1 on one of the terms, 


E E <^ F i ^ el - E fj i<> = E E i^ eI V? > 

i ci j i a v — 


=0 

Brillouiris theorem 






( 2 . 12 . 12 ) 


y 


=0 

orthogonal states 


so neither the single excitations contribute to the MP2 energy. Finally, the double excitations 
give the following term 


EE 

i ci 
j>ib>a 


<^ F i^ el -EAi<> 

1c 


E E QWeetyff) 

i a 


E E <XiXj\g\XaXb) - ( XiXj\g\XbXa > 

i ci 
j>ib>a 


(2.12.13) 
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FT 

±± 

FT 
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Figure 2.9: Illustration of the CASSCF method for a system of 8 electrons in 8 orbitals. The active space 
consists of 4 electrons in 4 orbitals, and in a CASSCF calculation we include all states where the 4 electrons 
are distributed over the 4 orbitals. According to eq. (2.12.5) this leads to A/qasscf = 8!/(4!) 2 = 70 in a CASSCF 
calculation whereas A/pci = 16!/(8!) 2 = 12870 for a full Cl calculation. 


so the MP2 energy becomes 



EL 

i a 


(XiXj\g\XaXb) ~ ( XiXj\g\XbXa > 
(ei - e a ) + (ej - e b ) 


2 


(2.12.14) 


where we in the denominator have used that the excitation energies for the unperturbed 
wavefunction is simply given as orbital energy differences. This can be simplified for a 
closed-shell system along the route taken to obtain eq. (2.6.39). 


2.12.3 Multiconfigurational SCF 

We introduce multiconfigurational SCF (MCSCF) methods by describing the complete active 
space SCF (CASSCF) method (Roos 1987). Based on a Hartree-Fock calculation, we select an 
active space and within this active space we allow for all excitations from the occupied states 
(see figure 2.9 for an example). As for the Cl method, the coefficients C A , in eq. (2.12.3) are 
determined variationally but in contrast to the Cl method, also the orbital coefficients c ; y 
in eq. (2.6.1) are reoptimized so the molecular orbitals are actually modified in an MCSCF 
method as compared to the Hartree-Fock method. MCSCF methods are thus suitable for 
molecules where the Hartree-Fock approximation is poor, for example for molecules with a 
nearly degenerate ground state or for describing the breaking of chemical bonds. However, a 
CASSCF calculation is by no means a black-box calculation. The selection of the active space 
requires a detailed insight about the molecular system and the cost of the calculation scales 
factorially with the size of the active space. 


2.13 Density functional theory 

Through the Hohenberg-Kohn theorem (Hohenberg and Kohn 1964), density functional 
theory (DFT) provides an alternative foundation to electronic structure theory as compared 
to the Schrodinger equation. DFT gives also a theoretical foundation to concepts of the 
electronic structure of molecular systems as electronegativity and chemical hardness. In 


Download free eBooks at bookboon.com 


58 














ATOMISTIC MODELS 


MOLECULAR QUANTUM MECHANICS 


addition, DFT based on the Kohn-Sham approach (KS-DFT) (Kohn and Sham 1965) has, 
essentially since the work by Becke (Becke 1993), become the main tool in computational 
quantum chemistry. The reason for its success is that KS-DFT includes electron correlation, 
and for many properties KS-DFT gives experimental accuracy, with a computational cost 
comparable to a Hartree-Fock calculation. 

The Hohenberg-Kohn existence theorem states that the ground-state energy and all other 
ground-state properties are uniquely defined by the electron density. That, for example, 
implies that the Hamiltonian in principle can be derived from the electron density. 
Mathematically this is expressed as an energy functional 3 of the electron density, E[p(f)], 

E[p(r)] = J V ext (f)p(f)dr+F[p(r )], (2.13.1) 

where V ext (r) is the external electrostatic potential where the most important term normally 
arises from the nuclei of the molecular system, and F[p(r)] is a so far unknown energy 
functional containing the remaining energy terms as the kinetic energy of the electrons and 
the electron-electron interactions. 


The number of electrons, n, is conserved, 


/ 


p(r)dr = n , 


(2.13.2) 


(2.13.3) 


and by applying the variational principle in section 2.7, the energy functional, E[p{r)], 
is minimized with respect to the electron density with the constraint that the number of 
electrons is kept constant as 

j 2 _(£[p(f)]-p/p(f)dr)= 0 , 

where 8 denotes a. functional derivative and p is a Lagrangian multiplier for the constraint in 
eq. (2.13.2). Eq. (2.13.3) is rewritten as 

' mpir)] ] 

< spin J, 


= E 


(2.13.4) 


text 


which may be regarded as the DFT equivalent to the Schrodinger equation. We put 14xt 
as a subscript to denote that it is kept constant (c.f. how we write partial derivatives 
in thermodynamics), e.g. the nuclear positions are kept fixed in the clamped-nucleus 
approximation as described in section 2.4. 


2.13.1 Electronegativity 


We identify p in eq. (2.13.4) as the chemical potential for electrons, and it also provides a 
definition for the electronegativity £ of the entire system (with a minus sign), 


■£ = /* = 



(2.13.5) 


3 A simple example of a functional: g\f(x)} - / fix) dx, i.e. we provide a function as the argument and the 

-OO 

functional returns a value, e.g. an energy. 
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The electronegativity is a concept used for a long time (Jensen 1996) and extensively in for 
example organic chemistry based on e.g. the electronegativity scales by Pauling (Pauling 
1932) and Mulliken (Mulliken 1935) to discuss reactivity in molecules, but it is first with 
DFT that the electronegativity obtained a proper definition (Parr and Yang 1989). Similarly, 
DFT provides the foundation for other concepts like the chemical hardness, Fukui indexes 
for reactivity in molecules, etc. (Parr and Yang 1989, Pearson 1997). We also use these 
concepts in the construction of a force field model based on electronegativity equalization 
in section 3.3.l.i. 


2.13.2 Kohn-Sham approach 

DFT has become the workhorse in computational quantum chemistry, and it is the Kohn- 
Sham approach that is used in the actual DFT calculations. In KS-DFT (Kohn and Sham 
1965), the energy functional F[p(r)] in eq. (2.13.1) is partitioned as 

F[p(f)] = J Vext(np(f)dr+FKE[p(f)] + F H [p(r)] + E xc [p(r )], (2.13.6) 

where FkeIp(c)] is the kinetic energy for an electron gas (ideal gas) but with the correct 
density, Fn[p(r)] is the Hartree term which corresponds to classical electrostatics and the 
Coulomb term in the Hartree-Fock approximation in section 2.10, and, finally, ExcAp (?)] 
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is the exchange-correlation functional and contains the remaining terms. We need a 
mathematical representation of the electron density, and the choice in the Kohn-Sham 
approach is to use orbitals, (p{?i), 


n 


pm = £ \<p(ji )\ 2 . 

i =1 


(2.13.7) 


We have thus left a 3-variable representation, as in p(r), and returned to a 3??-variable 
representation as in the molecular orbital approach. Furthermore, we can now use the same 
basis sets (see section (2.11.2) as in molecular orbital theory. Following eq. (2.13.4) for each 
term in eq. (2.13.6), we have for example, 


Vxc = 


' SExclpm] 

k 6p{r) 


text 


(2.13.8) 


where Vxc is the exchange-correlation functional. In the Kohn-Sham approach, eq. (2.13.4) 
leads to 

(hfCE + Text + Vr + hxc) <fii = £i<Pi ■ (2.13.9) 

Recall the Fock operator in eq. (2.10.9) 

fi = h i + it[Jj-'#])> (2.13.10) 

i =i 


where h, is the one-electron term here corresponding to Vre and V ext , Jfj is the Coulomb 
term which here corresponds to Vu, and the exchange term jfrj is replaced by the exchange- 
correlation functional Vxc- Thus it is reasonable to regard the terms in the Kohn-Sham 
approach in eq. (2.13.9) as a modified Fock operator, 

fP FT = VkE + Text + Vu + Vxc , (2.13.11) 


which includes an ad hoc correction for electron correlation in the exchange-correlation 
functional Vxc- The exact form of Vxc is unknown, but many functionals (and with a 
functional, we mean the choice of Vxc) have been suggested. Many of these functionals have 
been developed to describe certain properties and many functionals even contain empirical 
parameters fitted to experiments, so the choice of a functional in an actual calculation 
have to be based on a thorough literature study on DFT calculations on similar molecular 
systems. The success of computational DFT is, however, indisputable. With the proper 
choice of functional and basis set, results can be obtained that rivals experiments and with 
a computational cost that is considerably less expensive than other methods to include 
electron correlation. 


Exercises 

Ex. 2.1 — Show that the eigenvalues are real and that the eigenfunctions are orthogonal for 
a Hermitian operator. 
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Ex. 2.2 — Show that the normalization factor for an n-electron Slater determinant is 1 / \fn\. 

Ex. 2.3 — The energy of the hydrogen molecule in the molecular-orbital model (LCAO 
approximation) is regarded, (a) Give an expression for the energy of the ground state, (b) 
Give an expression for the triplet state with the lowest energy. 

Ex. 2.4— Show the relations for the antisymmetrizing operator si in eqs. (2.6.25) 
and (2.6.26). 

Ex. 2.5— Given the exact solution, J6\pi = Eip/j, show that for a trial function, xpi, the 
variation theorem is fulfilled, i.e., Eq > E$ where the equal sign holds when the trial function 
is identical to the exact solution. 

Ex. 2.6— Show the Rayleigh-Ritz method, i.e., apply the variational theorem on a trial 
function, 

Vo = E 

i 

expanded in a set of functions (pi . (a) Show first that the Rayleigh-Ritz method leads to the 
secular equation 

T. g | a it—E qShc^ = o 

i 

where = {(pi\^\(pk) and = {(pi\(pk)- (b) Rewrite this equation in terms of a secular 
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determinant 

|H-£ 0 S| = 0 

where is a matrix element of H, and similarly, S// v - is a matrix element of S. 

Ex. 2.7 — What is the dependence of the choice of origin on the molecular dipole moment 
for 

a) a neutral molecule ? 

b) a charged molecule ? 

Ex. 2.8 — What is the dependence of the choice of origin on the molecular second moment 
for 

a) an uncharged and unpolar molecule ? 

b) an uncharged and polar molecule ? 

c) a charged molecule ? 


Ex. 2.9— Showeq. (2.9.40) by using integration by parts (partial integration). 

Ex. 2.10— Show the Roothaan-Hall equations in eq. (2.11.8). 

Ex. 2.11 — Show the Gaussian product rule in eqs. (2.11.15) and (2.11.16). 

Ex. 2.12— Show that the molecular Hamiltonian in eq. (2.4.2) commutes with the Fock 
operator in eq. (2.10.9). 


Ex. 2.13 — 


Ex. 2.14 — 


Ex. 2.15 — 


Show eq. (2.8.37). 

Show eq. (2.8.38). 

Solve the Schrodinger equation 

d 2 xfj 
dcp 2 


in eq. (2.A.40), 
= -m 2 y/, 


and normalize it. 


2.A Quantum-mechanical model systems 

We here present some model systems in quantum mechanics which can be solved exactly. 
The selection is limited to model systems used in subsequent chapters. 

2.A.1 Translation - Particle in a one-dimensional box 

The model system is described in figure 2.10. It is a one-dimensional system where we in 
one region, 0 <x< a, have a potential energy that is zero, V = 0. In the remaining part, x < 0 


Download free eBooks at bookboon.com 


63 



ATOMISTIC MODELS 


MOLECULAR QUANTUM MECHANICS 



Figure 2.10: The one-dimensional box. 


and x> a, we have an infinite potential energy, V = oo, 


V{x) = 


0 0 <x< a 

oo x < 0 or x > a 


(2.A.1) 


The Schrodinger equation in one dimension is 


-h 2 d 2 y/ 
2 m dx 2 


+ Yy/ = Ey/. 


(2.A.2) 


In the region outside the box, x < 0 and x > a, the potential energy is infinitely high and 
it is therefore not possible for the particle to be in this region. The probability to find the 
particle is therefore zero in this region, y/ 2 = 0, and thus also the wavefunction is zero, 
y/ = 0. It gives two boundary conditions, y/{0) = 0 and y/(a) = 0, that is used to solve the 
Schrodinger equation for the central region. For the central region, 0 < x < a, where V = 0, 
the Schrodinger equation becomes, 


d 2 y/ -2mE 
dx 2 h 2 ^ 

This differential equation has a solution of the form, 


(2.A.3) 


y/{x) = Asin(ax) +Bcos(bx). 


(2.A.4) 


If the trial solution in eq. (2.A.4) is put into the Schrodinger equation in eq. (2.A.3), 

-2mE 


- Aa sin(ax) -Bb cos (fix) = 


h 2 


(Asin(ax)+ Bcos(fix)) , 


(2.A.5) 


it leads to 


a 2 = b 2 


2mE 

h 2 


(2.A.6) 
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Furthermore, using the boundary condition, = 0, and that the wavefunction is 

continuous leads to that B = 0 since cos(0) ^ 0. So far, we have achieved, 


if/(x) = A sin 


V2mE 
h 


-x . 


Using the second boundary condition, = 0, gives 


A sin 


V2 mE 
h 


a = 0, 


and using 
so that 

gives 


sin(w) = 0 

s/lmE 

h 


u = ±nn, n = 0, l,2,...,oo, 


a = ±nn, n = 0,1,2, ...,oo, 


mi 


a 


y/{x) = Asin +—x = +Asin —x , n = 0,l,2,...,oo, 


nn 


a 


(2.A.7) 


(2.A.8) 

(2A.9) 

(2.A.10) 

(2.A.11) 


since sin(-x) = - sin(x) .A is a constant with arbitrary sign, so we can simply replace + A with 
A. For n = 0, we note that y/(x) = 0. This solution is not allowed since the probability to find 
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Figure 2.11: The wavefunctions, (//,, (x), (in blue) for the one-dimensional box 


the particle, i/r 2 (x), becomes zero in the entire space, and the particle has to be somewhere. 
For the region, 0 < x < a, the wavefunction thus becomes 

(nn\ 

t/r(x) = AsinI—xl , n = 1,2,...,oo . (2.A.12) 

Normally, we assume that the wavefunction is normalized, 


J if/ 2 (x) dx= 1 
o 


[2 (nn . 

t/r(x) = y- sin I —x | , n = 1,2,..., oo . 


(2.A.13) 


Using eq. (2.A.10), the energies are obtained as 

n 2 n 2 h 2 n 2 h 2 

En — 


2ma 2 8ma 2 ’ 
The excitation energies, A E n , are thus given as 


n = l,2,...,oo. 


(2.A.14) 


A E n = E n+ 1 - E n = (2 n + 1) 


h 2 


8 ma 2 ' 


(2.A.15) 


The wavefunctions for a few n are shown in figure 2.11 and the corresponding probabilities 
are shown in figure 2 . 12 . 


2.A.1.1 Particle in a two- and three-dimensional box 

The model system for a two-dimensional box (can be used as a model for a surface) is a trivial 
extension of the one-dimensional box in section 2.A.I. The box has two side lengths, a and 
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Figure 2.12: The probabilities, ty 2 n {x), (in blue) for the one-dimensional box 


b } and the potential energy, V, is given as 


V (x, y) 


0 0 < x < a and 0 < y<b 

oo elsewhere 


(2.A.16) 


As for the one-dimensional box, the wavefunction, t//(x, y) = 0, is zero in the region where 
V = oo. For the region within the box, 0 < x < a and 0 < y < b, the Schrodinger equation is 


d 2 y/ d 2 y/ -2mE 
dx 2 + dy 2 H 2 ^ 


(2.A.17) 


As in i eq.2.1.5, we can use separation of variables and the wavefunction is given as the 
product of two one-dimensional functions, 


t \f/{x,y)=y/ x {x)y/ y {y) 


(2.A.18) 


Using the solutions of the one-dimensional box in eqs. (2.A.13) and (2.A.14), the wavefunc¬ 
tion becomes 


n x n 


\fab V ci 
and the energy is 


n y n 


i/r(x,y) = —j= sin| —^x |sin| ~^~y\ - n x = 1,2,...,oo , n y = 1,2,...oo, 


j n x n y 


8 m 


^ + n y 

a 2 


b 2 


A further generalization to three dimensions where 


V(x,y,z) 


0 0 < x < a and 0 < y < b and 0 < z < c 

oo elsewhere 


(2.A.19) 


(2.A.20) 


(2.A.21) 
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gives the wavefunction 

y/(x,y,z) = 


8 f n x n 

——sin -x|sm 

abc \ a 


n y n 


n z n 
y I sin I- z I , 


(2.A.22) 


which has the following possible quantum numbers: n x = l,2,...,oo ,n y = l,2,...oo , and 
n z = 1,2,... oo. The energy becomes 


h n x n y n z ~ ^ 


n z n v n‘ 

[he _ y_ £Zl 

a 2 b 2 c 2 


(2.A.23) 


2.A.2 Vibrations - Harmonic oscillator 


The harmonic oscillator may be regarded as a simple model for the vibration of a diatomic 
molecule around the equilibrium bond length, R e , where q = R - R e is thus the deviation 
from R e . The Hamiltonian is given as 




-H 2 d 
2p dq 2 


+ ^kq 2 , 


(2.A.24) 
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Figure 2.13: Particle on a ring with a radius r. The particle has the momentum p and the angular momentum 1. 


where /.< is the reduced mass for a diatomic molecule, 


mi m2 

h = -, 

mi + m2 


(2.A.25) 


and k is the force constant obtained as the second derivative of V{q) calculated in q = 0. In 
eq. (2.A.25), mj is the mass of atom I. The solutions for the Schrodinger equation are known 


1 f/ n {q) = A n H n (y)e~ 


(2.A.26) 


where e 2 is a Gaussian function, H n (y) are the Hermite polynomials, and 


y = 


a = 


(h 2 Y 

( 1 n \ 

— , 

A n = \an 2 2 n\\ 

w 


(2.A.27) 


The energies, e n , are given as 


1 


£ n = \ n + - \hv = \ n + - \flLo = \ n + 


1 


11 he 


2 A 


(2.A.28) 


which is expressed either with the frequency, v, the angular frequency, u>, or the wavelength, 

A. 


2.A.3 Rotation - Particle on a ring 

First we regard classical mechanics for a particle moving on a ring with a radius r = |r|, see 
figure 2.13. The particle has a momentum p = mv, where m is the mass of the particle and v 
is its velocity, and we define the angular momentum 1 as 

1 = r x p . (2.A.29) 


The velocity v = |v| can be written as 


v = 2nrv = rco , 


(2.A.30) 
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where v is the frequency (number of rotations per time unit) and co = 2nv is the angular 
velocity (radians per time unit). For l = |1| and p = |p|, we get from eq. (2.A.29) that l = rp. 
Since p = mv = mrco, we get l = mr 2 a> = Ia> where we define the moment of inertia I of the 
particle as 

/ = mr 2 . (2.A.31) 

The kinetic energy K, as well as the total energy E since the potential energy is 0, can thus be 
written as 

ll / 2 

E = K= -mv 2 = -Ico 2 = — . (2.A.32) 

2 2 2 / 

We will now demonstrate that the angular momentum 1 is quantized. A particle can rotate 
both ways, ±p, thus 1 in figure 2.13 can be directed both upwards and downwards, 

l z = +pr . (2.A.33) 


According to de Broglie each particle behave like a wave as 

A- h -, 

P 

where A is the wavelength. Thus the angular momentum becomes 

hr 


l z = ± 


A 


(2.A.34) 


(2.A.35) 


If we rotate the angle 2 n, we need to get the same l z , thus the number of wavelengths per full 
rotation have to be an integer mi, i.e. ra/A = 2nr, leading to 


hr 

h = +y =hmi, 


mi = 0 ,± 1 ,± 2 ,... , 


(2.A.36) 


where as before h = h/2n and the + is included in m/. Note that mi = 0 corresponds to an 
infinitely large wavelength A. The energy of the system in eq. (2.A.32) becomes 

' 2 

(2.A.37) 


/2 m 2 ,h 2 

E= — = —— 
21 21 


where we have degeneracy for all m/ apart from m/ = 0. The Schrodinger equation for the 
motion in the xy-plane (no z component) is given as 

- h2 (^ + Jt) V/ = 


jt// = Ey/ 


(2.A.38) 

for the case when the potential energy is zero. For a particle on a ring, we use cylindrical 
coordinate, r and (p , where the radius r is fixed (see figure 2.14), leading to the Schrodinger 
equation, 

-h 2 1 d 2 y/ 


= Ey/. 


(2.A.39) 


2m r 2 dcp 2 

Using eq. (2.A.37) for the energy and the definition of the moment of inertia in eq. (2.A.31) 
leads to 

d 2 y/ ? 

- 7 r = -m,y/. (2.A.40) 


dcp 2 


Thus the wavefunction becomes 


= 


3 i mnp 


\Z2n 


(2.A.41) 


which is normalized (see exercise 2.15). 
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Figure 2.14: Coordinate system for cylindrical coordinates 


2.A.4 Particle on a sphere 

For a particle on a sphere, we adopt spherical polar coordinates, r, 0 and (p (Arfken et al. 
2013, ch. 3), 

x = rsin(0)cos(<p) 0 <d <n 
y = r sin(0) sin(<p) 0<(p<2n 

z = r cos (0) r is a constant 


As for a particle on a ring, the potential energy is zero resulting in the following Schrodinger 
equation, 

_ 

-V 2 y/{r,e,(p) = Ey/{r,6,(p) 


2m 


For spherical polar coordinates, 


where A 2 is the Legendrian, 


o 1 d 2 1 2 

V :? - = - — r + ^A ‘ 2 , 
r dr z r z 


Id 2 Id 


(2.A.42) 


(2.A.43) 


+ 


sin 2 [0)d(p 2 ' sin(0) 66 50 ' 

The Legendrian has known eigenfunctions and eigenvalues a 


(2.A.44) 


A 2 Yi m (0, cp) = -/(/ + 1) Y lmi (0, cp ), (2.A.45) 

where Yi mi {9,(p) is a spherical harmonics which are tabulated (see e.g. (Arfken et al. 2013, 
ch. 15)). The eigenvalues / and m/ are quantized as, 


• / - 0,1,2,...,oo 

• ini = 0 , ± 1 ,... ± / . 
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Note that the allowed values of m\ depends on /, which is the reason for its notation (used 
also in the previous section). For a particle on a sphere, r is constant, so we get for the 
Schrodinger equation, 

~2j ^ Imi = , (2.A.46) 

so that 

E lmi = lil + l)jj. (2.A.47) 

The energy Ei m depends only on l so we have a degeneracy factor g = 21 + 1, i.e. the number 
of degenerate states. 


2.A.5 One-electron atom 


We will here solve the Schrodinger equation for the one-electron atom discussed in 
section 2.5.1 where atomic orbitals were introduced. Within the Born-Oppenheimer 
approximation (see section 2.4) we regard a nucleus at rest in a fixed position. The potential 
energy operator Y is given from eq. (2.5.1) as the Coulomb interaction between the nucleus 
and the electron, 


f = 


-Ze 2 

4neor 


(2.A.48) 
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Since Y only depends on one variable, the distance between the nucleus and the electron 
r, it is tempting to use spherical polar coordinates ( r,0,(p ) as in appendix 2.A.4 instead of 
Cartesian coordinates (x, y, z ). The Schrodinger equation becomes 


-h 2 

- + Yw = Ew 

2m e r r r 


(2.A.49) 


where m e is the mass of the electron and V 2 in spherical polar coordinates was given in 
eq. (2.A.43). The wavefunction y/{r,6,q)) is separable as 


i rir,o,(p) = Rirme,<p). 


(2.A.50) 


Thus the Schrodinger equation in spherical polar coordinates becomes 

^-(-§ I [rR{r)Y{e,(p))] + \A 2 R{r)Y{e,(p)-^-R{r)Y{e,(p) = ERir)Y{e,(p). (2.A.51) 

The angular part is identical to the solution for a particle on a sphere in eq. (2.A.45), 

A 2 Yi m (0, q,) =-1(1 + 1) Y lmi (0, <p) , (2.A.52) 

where Yi mi {6,cp) were the spherical harmonics. For the radial part, we get the following 
equation 

-h 2 (1 ,2 , A 1 . , , -Ze 2 

-— \-h rR{r ) + -x -/(/4-1) R(r ) - -- R(r) = ER(r ), (2.A.53) 

2 m e \r 0r )r z 4ji£ 0 r 

for which the solutions are the associated Laguerre functions, R n i(r) (see e.g. (Arfken et al. 
2013, ch. 18)). The wavefunction thus becomes 


y/nlmi Y ) — Rnl (^) T imi (0> ty) > (2.A.54) 

which is the result used in eq. (2.5.2). 
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3.1 Introduction to force fields 


In a force field, the potential energy surface of a molecule or a system of molecules is 
described with a relatively simple model based on atom-type parameters and approximate 
analytical expressions for a set of energy terms where each energy term has a physical 
interpretation. The motivation for constructing a force field is in particular that we in 
molecular dynamics (MD) simulations need to calculate the forces between a large amount 
of particles (perhaps 10.000 atoms) repeatedly. A typical time-step in an MD simulation is 1 fs 
so for an MD simulation of 100 ns we need to calculate all the interatomic forces 10 8 times. 
A standard force field for the potential energy surface, V(R 3N ), may look like 


bonds 


angles ^ 


torsions 


v '( r3N ) = TTf ('.-M 2 +Xf (e,-M 2 +"lTf (i + cos(w-ri)) 


bond stretches 


angle bends 


torsional motion 


+ Jl 4£ ij 

i,j>i 


(l„ y 

I" 

\6\ 

| 


1 " 



( li ( lj 
471£qR 


(3.1.1) 


ij 


Lennard-Jones term 


Coulomb term 


intermolecular interactions 


where each energy term will be described in some detail below. Each energy term has 
a physical meaning, but is in this chapter introduced in a phenomenological way based 
on empirical experience. Each term contains variables related to the molecular geometry, 
e.g. the bond length the bond angle, 0/, the torsional angle, ojj, and the interatomic 
distance, Rj j. In addition, each energy term contains atom-type parameters, e.g. the force 
constant, ki, and the equilibrium bond length, Z ; -,o, for the bond stretching term. Even if it 
is desirable to have only one value of each atom-type parameter for each element, it is in 
many cases not possible as for example for the carbon atom since a carbon atom in a methyl 
group, -CH 3 , and in a carbonyl group, >C=0, has quite different properties (e.g. the carbon 
atomic charges, qi, have different sign). The transferability of atom-type parameters, i.e. the 
applicability of the same set of values of the atom-type parameters to a variety of different 
types of molecular systems (e.g. proteins, nonpolar polymers, ionic crystals, etc.) is thus a 
central issue when discussing force fields. 

In the construction of a force field, we have to make a choice of which energy terms to 
include, the functional form of each energy term and how to obtain the values of the atom- 
type parameters. Historically, force fields have been parametrized against experimental 
data, termed empirical force fields, as for example structural and thermodynamic data. 
Semi-empirical force fields are still common, where usually the atomic charges in the 
Coulomb term in eq. (3.1.1) are obtianed from quantum chemical calculations and the 
remaining energy terms are parametrized from experimental data. Finally, a force field 
can be parametrized from a set of quantum chemical calculations on model systems. 
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Quantum chemical data is in principle preferred since it provides a consistent data set on 
the microscopic level, whereas most available experimental data are on the macroscopic 
level, i.e. measured at a given temperature and pressure. The limitation of the data set used 
in the parametrization limits the applicability of the force field, e.g. if only data at room 
temperature and ambient pressure are included in the parametrization it cannot be expected 
to be applicable at other thermodynamic states. 

We will in this chapter give an introduction to various force field terms based on a 
phenomenological approach, i.e. the energy terms are introduced in ad hoc manor and not 
derived from an underying and more fundamental theory. However, in section 3.4 some of 
the contributions will be derived from quantum mechanics. We divide the force field into 
two types of terms, bonded and non-bonded interactions. With bonded interactions, we 
refer to terms that describe covalent bonding, whereas non-bonded interactions describe 
relatively weak through-space interactions. Non-bonded interactions are not restricted 
to intermolecular interactions, but also includes interactions in relatively large molecules 
between atoms in different parts of the molecule. For example, atom 1 in the larger of the 
two molecules in figure 3.1 cannot see the difference between the interaction with atom N 
in the same molecule and the interaction with the atoms in the small molecule, so these 
interactions need to be treated on an equal footing in a force field model. 
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4 



N- 1 


Figure 3.1: Sketch to illustrate that the interactions 
between atoms of different parts of the molecule, 
in this case atom 1 and atom N, are essentially the 
same as the interaction between atom 1 and the 
atoms in the small molecule. 



Figure 3.2: The Morse potential for a diatomic 
molecule. D e is the dissociation energy and R e is the 
bond distance at the energy minimum, respectively. 


3.2 Force-field terms for covalent bonding 

We will here introduce the most common force-field terms for describing covalent bonding: 
bond stretching, angle bending and torsional motion. 

3.2.1 Bond stretching 

Let us regard the potential energy, V{R), as a function of the bond distance, R, of a 
diatomic molecule which often can be represented accurately with a Morse potential (see 
figure 3.2) (Morse 1929), 

V(R) = D e (< s ~ 2a ( R - R e ) - ) , (3.2.1) 

where R e is the equilibrium distance, i.e. the bond distance at the minimum of the potential 
well, D e is the depth of the potential surface (dissociation energy), i.e. minus the potential 
energy at R e , and a is a parameter describing the width of the potential well. The Morse 
potential can be written in different ways and is here given so that the potential energy 
approaches zero at infinite separation. 

As a first approximation, we represent a covalent bond with a classical harmonic oscillator, 
i.e. two atoms are connected with a spring with a force constant, fc.We make a Taylor 
expansion of the potential energy, V{R), around the equilibrium geometry, R e . If we 
introduce the Dunham expansion parameter Q = {R- R e )/R e (Dunham 1932a;b), the Taylor 
expansion around R = R e , i.e. Q = 0 becomes 

V(Q)= V(0) + V' (0)Q + -V"(0)Q 2 +-F (3) (0)Q 3 + — F (4) (0)Q 4 ... , (3.2.2) 

zero level V r '(0)=0 in the minimum ^ ^ v ' s “ v ' 

Harmonic term Anharmonicity Quartic term 

where the zero level, V (0), can be ignored (the potential energy surface can always be shifted 
with a constant energy). The linear term in Q is 0 since the gradient is 0 at the minimum. The 


Download free eBooks at bookboon.com 


78 








ATOMISTIC MODELS 


FORCE FIELDS 



Figure 3.3: A diatomic molecule represen¬ 
ted as two atoms connected with a spring. 
According to Hooke’s law, the force as a 
function of the bond distance, F{R), is given 
as F{R) = -k{R- R e ), where lc is the force 
constant and R e is the equilibrium bond 
length. 



Figure 3.4: A Taylor expansion of the Morse 
potential. The full line denotes the Morse 
potential, whereas the dash-dotted (red) 
line shows a truncation after the harmonic 
term, the dashed line (green) shows a trun¬ 
cation after the anharmonic (cubic) term, 
and the dotted line (blue) shows a trunca¬ 
tion after the quartic term, respectively. 


quadratic term in Q is the harmonic term, we denote the cubic term in Q the anharmonicity 
of the potential energy, and we finally have the quartic term. If the anharmonicity and 
other higher order terms in the Taylor expansion are ignored, we get the classical harmonic 
oscillator, 

lc 

V(Q) = -Q 2 with k = V" (0) . (3.2.3) 

A Taylor expansion of the Morse potential in eq. (3.2.1) is shown in figure 3.4. Truncation 
after the harmonic term gives a reasonable description around the minimum. At large 
separation between the atoms, Q — oo, also V(Q) -*■ oo, so the dissociation limit is 
not described correctly. If the model, however, is restricted to describing non-reacting 
molecules where the bond length always is expected to be close to R e , truncation after 
the harmonic term is a reasonable approximation. Truncation after the anharmonic 
term (again see figure 3.4) leads to that the energy V(Q) — -oo for large Q. Whereas a 
harmonic approximation may be useful in the sense that it keeps the molecules in a stable 
configuration, a truncation after the anharmonic term may be disastrous in an energy 
minimization scheme or in a molecular dynamics simulation where the energy cannot 
diverge towards -oo. A truncation after the quartic term has essentially the same features as 
a truncation after the harmonic term but it resembles the Morse potential more accurately 
at slightly larger Q than the harmonic approximation. 
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The Simons-Parr-Finlan (SPF) expansion, where the expansion parameter Q is chosen as Q = 
(i? - R e )/R (Simons et al. 1973) is shown in figure xx and it has a superior performance with 
respect to convergence as compared to the Dunham expansion essentially since Q — 1 when 
R — oo for this choice of Q. A generalization is given by the Thakkar expansion (Thakkar 


1975), 


where 


V(R) = c 0 (p)A 2 


OO 

1+ £ C n (p) A" 

^ n -1 


HR, p) = sgn{p) 



{ 1 p> 0 

0 p = 0 

-1 p< 0 


(3.2.4) 


(3.2.5) 


In eq. (3.2.4), p is a nonzero real number. The Thakkar expansion reduces to the Dunham 
expansion for p = - 1, to the SPF expansion for p = 1, and to the Lennard-Jones potential 
that will be discussed in eq. (3.3.53) for p = 6 and c n {p ) = 0 for n > 1. 


3.2.2 Angle bending 

We use the water molecule as a typical example of an angle bending term (see figure 3.5), 
where the equilibrium angle, 6 e , between the O-H bonds is 104.5°. As for bond stretching, 
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Figure 3.5: Definition of the internal co- Figure 3.6: A double minimum potential for 
ordinates, the bond lengths, R\ and R 2 and the angle bending term, e.g. for the water 
the bond angle 6 , for the water molecule. molecule. 


we normally use a harmonic oscillator approximation for the angle bending term, 


^( 0 ) = y ( 6 >- 6> e ) 2 


(3.2.6) 


which is an approximation that is valid when the angle 6 is close to the equilibrium angle 


Q e . A better model is given in figure 3.6 for the water molecule. However, for example for 
the water molecule at 6 = n, which corresponds to a linear water molecule, the potential 
energy barrier is so high that the linear molecule will not exist unless we are at extremely 
high temperatures. It is therefore in most cases sufficient to adopt a Taylor expansion around 
one of the minima as in eq. (3.2.6). For ammonia, on the other hand, the umbrella motion 
that inverts the molecule cannot be ignored at ambient conditions and eq. (3.2.6) would be 
a severe approximation. 

3.2.3 Dihedral terms 

A dihedral angle is defined as the angle between two planes as illustrated in figure 3.7. If the 
atoms 1, 2 and 3 define one plane, and the atoms 2, 3 and 4 define the second plane, the 
dihedral angle co is the angle between these two planes. The dihedral (or torsional) energy is 
given as a periodic function 



(3.2.7) 


where n gives the periodicity (e.g. n - 2 for 180° periodicity and n - 3 for 120° periodiicy, 
respectively), V n is the rotational energy barrier and y defines where the dihedral angle is 0. 
The energy term can also be written as the real part of a Fourier series, 


V{a>) = y C m cosM” 


(3.2.8) 


n 


Some molecules, as for example benzene, are expected to be planar, and often we would like 
to add a constraint to keep the molecule planar. If we add constraints, for example by adding 
Lagrangian multipliers in a constrained energy minimization, the molecule is kept planar all 
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2 4 



Figure 3.7: Definition of dihedral angle. 


the time. Alternatively, we can add an extra energy term, referred to as an energy restraint or 
an improper torsion term, 

V(w) = k a) (l- cos2a>) (3.2.9) 

where the molecule is kept almost planar if k 0J has a very large value. 


3.2.4 Cross terms 


Sometimes cross terms are included describing the coupling between e.g. two bond 
stretches or a bond stretch and an angle bending term. Let us take the water molecule in 
figure 3.5 as an example where the intramolecular motion is described by two bond lengths 
R\ and R 2 and a bond angle 6. We carry out a Taylor expansion around R\ t o, i? 2 ,o and 6q, 


V[R l> R 2> e) = V[R l , Q} R 2 ,M + {R l -R l , Q ) 


dV 


dRi 


+ (i?2 - ^2,o) 


dV 




1 


Rl, 0 ,^ 2 , 0,00 
d 2 V 


dRo 


#1,O>#2,O>0O 


^1,0>^2,0><70 1 


1 , , 2 d 2 V 

+ 2(«2-R 2 .o) ^ 


+ (i?i - ( R2 - ^2,0) 


^1,0»^2,0»^0 
d 2 V 


+ 2a* 


^1,0^2,0>^0 

d 2 V 


Ri,0’R2,0’9q 


dR 1 dR 2 


+ {Ri-Ri,o) {6-60) 


d 2 V 


+ 


(i? 2 -i? 2 ,o) (6>-6>o) 


d 2 V 


3R 2 d6 


Ri,0’R2,0’9o 

+ ... 


dRidd 


Rl, 0^2, 0,00 

(3.2.10) 


Ri,0’R2,o>9o 


where the last three terms are coupling terms. It is common that only a few of the coupling 
terms are important. For water it is suitable to introduce a symmetric stretch coordinate, 


Qi - Qi,o - (Ri ~ Ri,o) + (R 2 ~ R 2 ,o) 


(3.2.11) 


and an antisymmetric stretch coordinate, 


Q 2 _ Q 2 ,o - (-Ri - Ri,o) - (R 2 _ ^ 2 , 0 ) 


(3.2.12) 
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where it turns out that only one of the cross terms, 


V{Qi,Q2 ,0) = ... + (Qi -Qi,o) [9-9q) 


d 2 V 

dQrfO 


+ ... 


Ql,O.Q2,O»0O 


(3.2.13) 


is of importance. When including cross terms, the atom-type parameters f?i,o, i? 2 ,o and 0 O 
(or Qi o, Q 2} o and 0 O ) do not correspond to the equilibrium geometry, which is realized by 
minimizing the potential energy surface in eq. (3.2.10). 


3.2.5 Summary of bonding terms 

We have from a phenomenological standpoint discussed the regular energy terms included 
in a force field to describe covalent bonds: bond stretching, angle terms and torsional 
motion, sometimes extended by coupling terms and additional terms for example describing 
hyperconjugation (how ^-conjugation affects bond stretching in neighbouring groups). 
Sometimes the following classification is used (Maple et al. 1994, Hwang et al. 1994, Allinger 
et al. 1996): 

• Class I: only harmonic terms, not any cross terms 

• Class II: anharmonic terms and cross terms 

• Class III: also additional terms, e.g. to describe hyperconjugation 


3.3 Intermolecular interactions 

In section 3.2, we discussed interactions mainly describing covalent bonds referred to as 
bonded interactions. In this section, we discuss intermolecular interactions, i.e. weak 
interactions as compared to covalent forces, as for example dispersion interactions in 
liquid argon, hydrogen bonding in liquid water and ion-ion interactions in electrolytes. 
Furthermore, the energy terms used to describe intermolecular interactions are also adopted 
for long-range interactions within the same molecule, as the energy terms discussed in the 
previous section are restricted to interactions between an atom and atoms up to three atoms 
away within the same molecule (see figure 3.1, i.e. the bonded terms include interactions 
between atom 1 and atoms 2, 3 and 4). 

As is in the previous section, we will introduce the energy terms for intermolecular 
interactions in a phenomenological way, whereas the connection to quantum mechanics 
is treated in section 3.4. The four most important intermolecular interaction energies are: 

• electrostatic energy: interactions between ions and/or polar systems using classical 
electrostatics. 

• induction energy: arising from that the electron density of a molecule is polarized by an 
electric field from the surrounding molecules resulting in induced electric moments (e.g. 
induced dipole moments). 
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• short-range repulsion energy: essentially arising from that the Pauli exclusion principle 
needs to be fulfilled. 

• dispersion energy: arising from that the motion of the electrons is correlated. 


3.3.1 Electrostatic interactions 


In quantum chemistry, the charge distribution of a molecule with N atoms is represented by 
a set of nuclear charges, {Zj, J = 1,2,..., AT}, and an electron density, p(r), so that 

J p{j) dr = n , (3.3.1) 

where n is the number of electrons in the molecule. In a force field, the most simple way to 

represent the molecular charge distribution is by a set of atomic charges, {qi,I = 1,2.A/), 

in conjunction with Coulomb’s law for the interaction between atomic charges, 


<7 iRj 


N N 

v=y y 

1=1 J=I+1 47T£oRlJ 


(3.3.2) 


where Rij is the distance between atoms I and /. The key problem in obtaining atomic 
charges is therefore how to partition the electron density into atomic contributions. An 
atomic charge therefore consists of the part of the electron density assigned to the atom plus 
the nuclear charge of the atom. Also, an atomic charge does not have a unique definition 
since it is not a measurable quantity, and in the literature we find literally hundreds of 
approaches to calculate atomic charges. 

It is reasonable to expect that the atomic charges of a molecule should reproduce the 
molecular electric moments discussed in section 2.9.1. It is always imposed that the 
molecule has the correct total charge, 


N 

q mo1 = < £ j qi. (3.3.3) 

7=1 

For a neutral molecule, q mo] = 0, even round-off errors giving a small molecular charge 
different from zero would give large errors in the electrostatic energy because of the long- 
range 1 IR distance dependence for charge-charge interactions. For small molecules, it is 
reasonable that the molecular dipole moment, 

N 

aC 1 = E Ri r i.« > (3.3.4) 

7=1 

and the molecular quadrupole moment (see eq. (2.9.32)), 

N ^3 J A 

Q™* = J1 C R 1^2 R I>a R I,P - 2 Rl 'T R I>r^apj > (3.3.5) 

are well described by the atomic charges qi. 
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If atomic charges are not sufficient to represent the molecular charge distribution, there 
are several ways to improve the model. A systematic way would be to add atomic dipole 
moments, /ij, a , and atomic second moments, Qi iU p, (or alternatively atomic quadrupole 
moments, Qi iU p) which would improve the corresponding molecular moments as 

N 

+ /*/,«, (3.3.6) 

7=1 


Qafi 1 = H <7 I R I,aRl,p + Vl,aRl,p + Rl,aVl,P + Q/,a/3 > (3.3.7) 

/=1 

where the latter can be converted to the quadrupole moment by eq. (2.9.32). 

Alternatively, extra charges placed outside the atoms, so-called virtual charges, may be 
added for describing lone-pairs, e.g. in water, or ^--systems, e.g. in benzene. Again, the 
representation of the electrostatics in a force field in terms of distributed charges, dipole 
moments, etc. is not unique and many models exist in the literature, but it appears like 
atomic charges extended by atomic dipole moments and sometimes atomic quadrupole 
moments is a more general and systematic approach. 

Atomic charges are often obtained from quantum chemical calculations using various parti¬ 
tion techniques of the electronic charge distribution as for example Mulliken charges (Mul- 
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liken 1955), Hirshfeld charges (Hirshfeld 1977), etc., which traditionally often has been com¬ 
bined with experimental data to obtain atom-type parameters for the other energy terms 
in a semi-empirical force field. For small molecules, atomic charges may be parametrized 
from quantum-chemically derived molecular dipole moments and quadrupole moments 
(eqs. (3.3.4) and (3.3.5)), but the dipole moment for larger molecules is not a very meaningful 
property since it is a sum of numerous important contributions where each contribution 
describes the local electrostatics around e.g. a binding site and needs to be modelled 
accurately. 

One way to address this is to parametrize the electrostatic potential on a set of lattice points 
{aj,j = l,...,m} around a molecule, 


£ qi 

(paj = L-r~> ( 3 - 3 - 8 ) 

against a quantum-chemically derived electrostatic potential. It is an attractive approach 
since the electrostatic potential is probed locally around all important interaction sites of 
the molecule. One potential problem is that the electrostatic potential calculated close to 
the molecule does not follow the classical behaviour for point charges in eq. (3.3.8) since at 
short distances, the quantum-chemically derived electrostatic potential is calculated within 
the molecular charge distribution and short-range damping terms are required. In general, 
parametrizations like this are tricky because of redundant data and over-fitting so expertise 
in multivariate data analysis and chemometrics is highly demanded. 


Example 3.1: The dipole moment of HCN. 


We can in principle get the correct dipole moment of HCN by using partial atomic 
charges in two different ways as illustrated in figure 3.8, i.e. either by a charge transfer 
between C and N or between H and C. With the choice of 5 = 0.5 e, both possibilities 
give approximately the correct dipole moment by using eq. 3.3.4, = -2.81 D for 

the “CN” case and p™ 1 = -2.61 D for the “CH” case, respectively, compared to the 
quantum chemical value of /if 110 ' = -2.70 D, but the two cases can, for this example, 
easily be distinguished by calculating the quadrupole moment. The quadrupole 
moments calculated by eq. (3.3.5), 0™ 01 , become -0.44 B for the “CN” case and 5.47 B 
for the “CH” case, respectively, compared to the quantum chemical value of 1.92 B. 

A reasonable distribution of the partial atomic charges for HCN obtained from fitting 
the electrostatic potential around the molecule is given in figure 3.9. These partial 
atomic charges give a dipole moment p'J 10 ' = -2.67 D and a quadrupole moment 
0™ 01 = 1.80 B again using eqs. (3.3.4) and (3.3.5), both in good agreement with the 
respective quantum chemical value. 

The quantum chemical calculations have been carried out with the NWChem 
package (Valiev et al. 2010) adopting the PBE functional (Perdew et al. 1996) and the 
cc-pVDZ basis set (Dunning Jr. 1989), resulting in the bond lengths Rcn = 1.170 A and 
Rch = 1.085 A, respectively, in addition to the quantum chemical values given above. 


Atomic partial charges are to a very small degree transferable, e.g. there is not one set of 
carbon atomic charges that can be used for all types of carbon atoms in e.g. methyl groups, 
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z 


Figure 3.8: Two possible distributions of partial atomic 
charges in HCN. 


Figure 3.9: Partial 

atomic charges of HCN 
obtained from fitting the 
electrostatic potential 

around the molecule. 


aromatic rings, carbonyle groups, etc. Let us consider a simple example, F 2 and HF. The 
atomic charges in F 2 , c/f = 0, because of symmetry reasons whereas HF is highly polar with a 
dipole moment of 1.8265 D (Muenter and Klemperer 1970) described by the partial charges 
qp = -qH = -0.41 e for the bond length 0.9168 A (Mann et al. 1961). Another example is to 
consider the carbon partial charges in benzene and formaldehyde. In benzene, qc ~ -0.05 e 
wheras qc ~ 0.40 e in formaldehyde, again demonstrating the lack of transferability of atomic 
partial charges. 


3.3.1.1 Electronegativity equalization model 


In the electronegativity equalization model (EEM), it is assumed that each atom in a 
molecule can be described by an atomic electronegativity, £/, and an atomic chemical 
hardness, qi (Mortier et al. 1985). If the electronegativities of two atoms are different, charge 
will flow from one atom the other until the molecular electronegativity (chemical potential) 
is the same everywhere (equalized). In addition, a work is required to charge an atom, which 
is determined by the chemical hardness (or capacitance) of the atom. Electronegativity and 
chemical hardness are central concepts in density functional theory and has been discussed 
thoroughly (Parr and Yang 1989, Geerlings et al. 2003). The molecular energy, V, for a 
molecule with N atoms becomes 


N i N 

v = + 


N 




hW 


q 


mol 



(3.3.9) 


where the third term is the regular Coulomb interaction between two atomic charges, qi 
and qj. The last term is a constraint that preserves the molecular charge, q mo1 , where /r 
is a Lagrangian multiplier. The atomic charges are obtained by minimizing the molecular 
energy in eq. (3.3.9) with respect to the atomic charges and the Lagrangian multiplier. We 
thus obtain 


dV 

dqp 


N 

= 0 = £k + VKqK + T^qp-q, 


L^K 


(3.3.10) 


and 


dV 

dq 


= 0 = q 


mol 



(3.3.11) 
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where the last term is simply the applied constraint. A set of N + 1 coupled linear equations 
is obtained which can be solved by standard numerical techniques 


m 
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(3.3.12) 


or alternatively the atomic charges are calculated as 
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(3.3.13) 


Thus, the atomic electronegativity, £/, and chemical hardness, rji, are atom-type parameters 
that has to be determined. The method has been parametrized in different ways, e.g. the 
Delft molecular mechanics (DMM) force field has been developed for hydrocarbons by 
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Donor Acceptor 



Figure 3.10: Water dimer. The angles 9 and co refer to the notation in eq. (3.3.15), where the dashed line 
indicates the direction of an oxygen lone-pair orbital. The distance Rh-o in eq. (3.3.15) refers to the distance 
between the H atom in the hydrogen donor molecule and the O atom in the hydrogen acceptor molecule. 


parametrizing experimental data (van Duin et al. 1994). Electronegativity equalization is 
useful in organic chemistry to determine where it is likely that electrophilic and nucleophilic 
attacks occur in a molecule since these types of reactions highly depend on the local 
differences in the electrostatic potential around the molecule. 


For a diatomic molecule, using the constraint qi + qj = 0, we get the solution (see 
exercise 3.1), 


-911 = 911 = 


jWi 

q I + qj-2Tf] ' 


(3.3.14) 


The sign of the partial charge is thus determined by the difference in electronegativity 
between the two atom-types in the dimer. Consequently, the model works for the example 
discussed above, F 2 and HF. The atomic partial charges qp in F 2 are zero by cancellation of 
~ in eq. (3.3.14), wheras they are non-zero in HF since £,h~^,f ^ 0. The magnitude 
of the atomic charges depends on all atom-type parameters as well as the bond distance in 
Tfj , so the EEM also includes a geometry dependence of the atomic charges important in 
e.g. modeling conformational energies in molecules. 


3.3.1.2 Hydrogen bonding 

Hydrogen bonds are in general difficult to model in a force field and extra terms are 
sometimes added to the force field as in the YETI force field (Vedani 1988), 


Vyeti - 


C 


pl2 p 10 

^H-O ^H-O 


cos 2 6 cos 4 CO 


(3.3.15) 


where A and C are parameters and the notation is explained for the water dimer in 
figure 3.10. This energy term is highly dependent on the angles 0 and co with the purpose 
to model that hydrogen bonds are strongly orientation dependent. 

Since hydrogen bonds are to a large extent electrostatic interactions (Stone et al. 1997), a 
more systematic way is to include atomic dipole moments in the model. A dipole-dipole 
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interaction can according to eq. (2.9.35) be written as 


*w=- 


3 (pj • Ru j [Ru ■ pj j - (pj ■ pj) 


R 


IJ 


(3.3.16) 


gives an appropriate dependence on the mutual molecular orientation as seen from the 
vector scalar products in eq. (3.3.16). 


3.3.2 Electronic polarization 

The electronic density of a molecule is polarized by an electric field from the electric 
moments of the surrounding molecules, normally referred to as the induction energy in a 
force field. We recall from section 2.9.2 that the polarizability ct a p is defined as the linear 
response to an electric field, Ep, 

^ d = a a pEp (3.3.17) 

where /i" ld is the induced dipole moment. We generalize eq. (3.3.17) to atomic polarizabilit- 
ies, oci^pj as 

= a riali E\* (3.3.18) 

where is the atomic induced dipole moment and E^ is the total electric field at atom 
I, i.e. it includes for example the electric field from an applied external field and the atomic 
charges of the surrounding molecules as well as the electric field from all surrounding atomic 
induced dipole moments. 

The discussion in section 2.9.2 for a classical polarizability in an external field is extended to 
a system of many polarizable particles. In addition to the electrostatic energy and self-energy 
given in eq. (2.9.41) for a single particle, a dipole-dipole interaction energy (see eq. (3.4.8)) is 
also included. The induction energy \4id of a system of polarizable particles in an external 
electric field, E is given as (Vesely 1977) 


fmd = - 


1 

2 



(3.3.19) 


where pp^ is the induced dipole moment of particle I, Tfj a p, is the interaction tensor 
of rank 2 defined in eq. (2.9.22), and V/ )Se if is the self-energy of particle I as defined in 
eq. (2.9.40), 

^,seif = \ («/, ap)' 1 R7 d a pTp ■ ( 3 - 3 - 2 °) 

The external electric field may arise from permanent electric moments of the surrounding 
particles of the system or other external sources. The system is polarized in such a way that 
the induction energy is minimized, 


dkind 


d/d 


ind 

K,j 



(3.3.21) 
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which gives 


Vk% = 


E K,y + 'll T KJ, r pL l j,p 
J*K 


N 


< 2 ) 


,ind 


(3.3.22) 


The atomic induced dipole moments are thus given from a set of 3 N coupled equations and 
is a true many-body term. If polarizabilities are included in a force field, the calculation of 
induced dipole moments will be relatively time-consuming as compared to other energy 
terms. The expression of the potential energy in eq. (3.3.19) may be simplified by using 
eq. (3.3.22) as 
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p ext 
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(3.3.23) 


The self-energy thus cancels the induced dipole-induced dipole interaction and half of 
the interaction between the induced dipoles and the external electric field. By inserting 
eq. (3.3.22), eq. (3.3.23) maybe rewritten as 
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(3.3.24) 
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Figure 3.11: Polarization in a system divided into two subsystems. 


Although the induction energy in eq. (3.3.19) is quadratic in the induced dipole moment, the 
final result in eq. (3.3.23) is only linear in . The result in eq. (3.3.23) is, however, only valid 
for the optimum p 1 "^ fulfilling the requirement in eq. (3.3.21), dl/ ind /dp“ = 0. Eq. (3.3.24) 
can therefore, for example, not be used to calculate the gradient of the induction energy. 


3.3.2.1 Distributed polarizabilities 

If a system is divided into subsystems (see figure 3.11 for an illustration), as for example a 
molecular polarizability into atomic contributions, there are two major contributions to the 
polarizability. We have a monopole contribution arising from charge transfter between the 
subsystems as a response to the difference in electrostatic potential between the subsystems, 
which we will describe by an extension of the electronegativity equalization method in 
section 3.3.1.l. Secondly, we have a local polarization within each subsystem where the 
leading term is described by an atomic induced dipole moment. 


3.3.2.2 Electronegativity equalization methods 


The electronegativity equalization method (EEM) for calculating atomic charges in sec¬ 
tion 3.3.l.i is extended with the interaction with an external potential, </>® xt . The molecular 
energy V in eq. (3.3.9) becomes thus 



Following the same procedure as in section 3.3.1.l, eq. (3.3.12) becomes 
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(3.3.25) 


(3.3.26) 
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Dividing q tot into the atomic charges defined in section 3.3.1.1 and induced atomic charges, 
g* nd , arising from the response to external electrostatic potential gives for g* nd , 


f n ind \ 

°l 1 
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T (0) 

1 12 
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••• 1 IN 
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-1 
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0 J 


or 


N 


qf d = Z Ai,<pf 

J= 1 


(3.3.28) 


where A\j is a relay matrix determining charge flow to atom I from an external electrostatic 
potential at atom /. 

If a homogeneous external electric field, £^ xt , is regarded, its electrostatic potential at atom 
/, <pf x is 

(pf- = -Rj.pEf' ■ (3.3.29) 


The response to E ^ xt thus becomes 


d qf A * 

g£s = -Z A ” R W- (3 ' 3 ' 30) 

We recall from eq. (2.9.36) that the molecular polarizability is defined as 




d^_ 

dEf x ’ 


(3.3.31) 


and since 


N 


fi nd =E . 

1=1 


we get 


N 

a a/5 = ~ Y R I,aMjRj,p ■ 
I,J=1 


(3.3.32) 


(3.3.33) 


For planar molecules, however, only an in-plane polarizability is obtained, and EEM can thus 
only give a part of the molecular polarizability. 


3.3.2.3 Point-dipole interaction model 

The point-dipole interaction model (PDI) model has many similarities with the electroneg¬ 
ativity equalization model in the previous section 3.3.2.2 in the sense that it relies on a set of 
native atom-type parameters that couple with each other through electrostatic interactions 
to give the molecular property. If we have a molecule of N atoms in an external electric field, 
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Ep^, and an atom-type polarizability, ap, 
moment, p“ ld , is given as 

,,ind _ 

M I,a ~ a I,af) 

where the last term is the electric field from all other atomic induced dipole moments in the 
molecule. Since we have 3 N coupled equations, the set of equations can be written in matrix 
form as 

/i ind = a (E ext + T (2) /i ind ) (3.3.35) 

which may be rearranged as 


is assigned to each atom, an atomic induced dipole 


N 


F^ext . 7^(2) in 

h I,P + L hj.apVj, 

m 


(3.3.34) 


l/ nd = (a -1 - T (2) ) 1 E ext = BE ext 
where B is the relay tensor defined as 

B= (a _1 -T®) -1 

The atomic induced dipole moment is given in terms of the relay tensor as 


N 

fii=LBiJ, a pE?p 

J= i 


(3.3.36) 


(3.3.37) 


(3.3.38) 
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It is noted that the molecular induced dipole moment is simply the sum of the atomic 


induced dipole moments, 


,ind . 


N 


/*« “ = E aC 

/=i 


(3.3.39) 


In addition, if a homogeneous electric field, Ej X ^ = £^ xt , is assumed, the molecular induced 
dipole moment is given as 


iC* = 


N 

B U,ctp 

1,1=1 


7ext 


(3.3.40) 


From its definition in eq. (2.9.36), 

|C* = (3-3.41) 

the molecular polariazability is obtained as 


N 

a a/3 1= E B U,a( J > 
I,J= 1 


(3.3.42) 


The definition of atomic polarizabilities is not unique, but one way to define it is as 


ind _ atom ^ext 
Pi,a ~ a I,ap n /3 


which results in 


a 


atom _ 
I,af} ~ 


N 

Jl B IJ,al5 

J 


Even if the atom-type parameter, ap, is isotropic, 


(3.3.43) 

(3.3.44) 


&P,a() ~ Mpdaf) 


(3.3.45) 


the resulting molecular polarizability includes the anisotropy because of the intramolecular 
electric fields. 

The PDI model thus gives a distributed representation of the molecular polarizability in 
terms of the relay tensor Bjj ta p. The relay tensor is interpreted as the response in terms of 
an atomic induced dipole moment at atom /, p 1 ^, arising from a perturbation by an electric 
field at atom /, Ej t p. The PDI model may also be regarded as a correction to an additive 
model by carrying out a Taylor expansion of the B tensor in eq. (3.3.33) around aT (2) = 0, 


B=(aT 1 -T (2) ) 1 a = |l + aT (2) + (aT (2) ) 2 + ...J a (3.3.46) 


which converges if aT (2) <1. It is also instructive to solve the two-atom problem, which 
results in (Silberstein 1917a;b) 


a \\ = 


at + ctj +4 (XiCCjlR 3 
1-4 oCiCijlR 6 


for the component parallel to the diatomic bond a y, and 


a± = 


at + aj - 2aiCLjlR 
1 - ajCtj/R 6 


(3.3.47) 


(3.3.48) 
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for the component perpendicular to the bond aj_. If R 6 -term is neglected, the isotropic part 
of the molecular polarizability is additive 

a mo1 = ^ (a,, + 2a±) = a; + ctj (3.3.49) 

which indicates that the isotropic part of the polarizability is easier to model than the 
individual tensor components and thereby also the anisotropy of the polarizability tensor. 


3.3.3 Dispersion and short-range repulsion 

If we consider the interactions in e.g. liquid argon, where electrostatic and polarization 
(in the absence of an external electric field) interactions are not present, the dominating 
interactions consist of attractive dispersion interactions and short-range repulsion. 


3.3.3.1 Dispersion 

Dispersion interactions arises from that the motion of the electrons are correlated, and 
London derived an equation for the dispersion energy from second-order perturbation 
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theory (see section 3.4.2 for details) resulting in (London 1930) 

Vdisp = ^ , (3.3.50) 

where Cq is a parameter with a postive value since the dispersion energy is always attractive. 
Higher-order terms in the perturbation expansion may be included, 


T/ _ ~C& Cs C io 

^disp - r6 r8 Rl0 + 

but normally only the R ~ 6 -term is included. 


(3.3.51) 


3.3.3.2 Repulsion 

The physical origin of the exchange-repulsion energy is the Pauli exclusion principle, i.e. 
that two electrons cannot be in the same quantum mechanical state. So when two closed- 
shell, e.g. two argon atoms, come close to each other they will repell each other just 
to fulfil the Pauli exclusion principle. The other energy terms commonly included to 
describe intermolecular forces (electrostatics, induction and dispersion), normally have 
quite significant approximations at short interatomic distances whereas the long-range 
behaviour is more or less correct. It is therefore common to define the repulsion energy 
Vrep more pragmatically as (Engkvist et al. 2000) 

Prep = Pint — Pele — Pjnd — Pdisp > (3.3.52) 

so that it contains also short-range errors in the other energy terms as well as higher-order 
energy terms that often are short-range, as e.g. for the dispersion energy in eq. (3.3.51). 
The repulsion energy Vr ep is therefore often parametrized from a set of quantum chemical 
calculations of the interaction energy Vj nt on molecular dimers and clusters and subtracting 
the other energy terms. 


3.3.3.3 Lennard-Jones potential 


Lennard-Jones suggested the following model for the interaction between two rare-gas 
atoms (Lennard-Jones 1931), 


/ 

Plj = 4e 

k 


a 

R 


12 



j4__Ce 
R 12 R 6 ’ 


(3.3.53) 


which consists of a short-range repulsion energy with an l? -12 -dependence and a more 
long-range attractive London dispersion energy as in eq. (3.3.50). We express the atom- 
type parameters either as e and cr or as A and Cq. The reason for the choice of the 
f? 12 -dependence is that the repulsive (a/R) ]2 term is the square of the attractive (cr/J?) 6 
term, which Lennard-Jones took advantage of in solving some statistical mechanical model 
systems analytically. 1 In a molecular force field, the Lennard-Jones potential is generalized 

1 The choice of the R ~ 12 -dependence is often attributed to computational efficiency. It is indeed 
computationally efficient, which may explain its present popularity, but Lennard-Jones presented his work 
in 1931 long before the first computers in the 1950’s. 
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to a pair-wise additive force field, 


N a N b 

= E E 4e u 


1=1j=i 



_ \ 6 ^ 
°u\ 


R u, 


(3.3.54) 


where Na and A/g are the number of atoms in molecules A and B, respectively. Eq. (3.3.54) 
contains atom-pair parameters, e// and ay/, which leads to M(M + l)/2 number of 
parameters for M atom types. It is favourable to reduce the number of parameters in the 
model to M parameters by applying the Lorentz-Berthelot mixing rules, 

a IJ =-[ a II + a Jj) ’ e IJ = V £ II £ JJ (3.3.55) 


Example 3.2: Interaction between an Ar atom and an ion. 


Is the leading interaction energy between an Ar atom and an ion a dispersion energy 
or an induction energy? It is an induction energy, and the argument is as follows. The 
ion with a charge q gives an electric field at the Ar atom, 

where we in the last step have assumed that both particles have been placed along 
one of the coordinate axis x, y or z. The Ar atom is polarized by the field and gets an 
induced dipole moment, 

Ba = a af)Ep 

We recall the induction energy in eq. (3.3.23) 

Mnd = — 2 Ra R<x = ~2^ a ^^^ a 

The electric field from a charge has an i? _2 -dependence so the induction energy has 
an i? -4 -dependence as compared to an R ~ 6 -dependence for the dispersion energy. 


3.3.3.4 Many-body interactions 

A pair-wise additive force field is of the form 

N N 

V =L E vij ( 3 - 3 - 56 ) 

i=i j=i+ i 

where Vjy is the pair-interaction energy. The electrostatic energy discussed in section 3.3.1 
is a true pair-wise additive energy term (see Coulomb’s law in eq. (3.3.2)). This pair-wise 
additivity is, however, lost if the atomic charges are calculated by the EEM in section 3.3.l.i 
since the values of the atomic charges depend on the surroundings. The induction energy 
in section 3.3.2 is a true A-body term in the sense that the induced dipole moment in 
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eq. (3.3.22) depends on the electric field of the induced dipole moments of all surrounding 
molecules leading to a set of coupled equations. 

The repulsion and dispersion energies are in the Lennard-Jones potential in section 3.3.3.3 
given as a pair-wise additive energy. The repulsion and dispersion energies are, however, 
only approximately pair-wise additive and corrections can be included as three-body, four- 
body, etc. energy terms, 

N N N N N N N N N 

V =L L %+E E E v UK+L E E E V IJKL + ... (3.3.57) 

1=1 J=I+1 I=1J=I+1K=J+1 I=1J=I+1K=J+1L=K+1 


A famous example is the three-body dispersion (Axilrod-Teller) term given here for a system 
of three identical particles as (Axilrod and Teller 1943) 


V AT = C 


3 cosyi cosj 2 cosy 3 + 1 

p3 d3 d 3 
ri 12 ri 23 ri 31 


(3.3.58) 


where the angles and distances are defined in figure 3.12 and C is a positive number 
proportional to Vjpa 3 where Vjp is the ionization potential and a is the polarizability. 
The Axilrod-Teller term, Vat, is for a system of rare-gas atoms around 8% of the London 
dispersion term in eq. (3.3.50) and should not be ignored in a study of e.g. liquid argon. 
For a system of polar molecules, e.g. liquid water, the London dispersion energy is only 
around 10% of the electrostatic energy so the Axilrod-Teller term can in this case to a good 
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Figure 3.12: Three-atom system with the notation for the Axilrod-Teller term. Rij in eq. (3.3.58) are the 
distances between the atoms. 

approximation be ignored since it contributes with less than 1% of the total interaction 
energy. 

It is a poor idea to approximate the induction energy in section 3.3.2 with three-body and 
possibly four-body terms as demonstrated in exercise 3.4, instead one should rely on the 
expression based on classical electrostatics discussed in section 3.3.2. 


3.3.4 Effective force fields 


With effective force fields, we mean force fields that effectively include some energy terms 
into some of the other energy terms simplifying the force field further. Effective force fields 
thus include the essence of the interactions needed for modeling a particular type of systems 
but their applicability area will be rather restricted. As discussed in eq. (3.3.56),, a pair-wise 
additive force field is of the form 

N N 

V = L L v u ( 3 - 3 - 59 ) 

/=i/=/+i 

where Vjj is the pair-interaction energy between atoms I and J. For N atoms, we thus have 
to compute N{N -1) /2 pair energies Vjj, and a major saving in computational time may be 
achieved from reducing the number of interaction sites N rather than simplifying Vjj. This 
is particularly true if the most expensive term in the calculation of Vj j is to calculate 


1 _ 1 

RlJ \J (*/ - -hO 2 + (y/ - y ,) 2 + (zj - Zj ) 2 


(3.3.60) 


where the divide and square-root operations are computationally expensive. 

An important example is electrolytes where it is common to replace the solvent molecules by 
a dielectric constant (relative permittivity), e, that shields the interactions between the ions 
in the calculation of the electrostatic energy, V e \ e , 


qiqj 


M M 

^ele = T T , D 
^nsoeRu 


(3.3.61) 


where M is the number of ions in the solution. Typical values for e is e ~ 2 for unpolar 
hydrocarbons and e « 78 for liquid water. 
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A common example is effective pair potentials where all many-body terms, including also 
electronic polarization discussed in section 3.3.2, are described by a pair-wise additive 
atomistic force field, 

V=V de + V LJ , (3.3.62) 

consisting of an electrostatic term Vei e in eq. ( 3 . 3 . 2 ) and a Lennard-Jones term Vy in 
eq. (3.3.54). For example for a pair potential of water, the atomic charges should for an 
effective pair potentential not reflect the gas-phase dipole moment of 1.85 D (Clough et al. 
1973, Dyke and Muenter 1973), but rather the dipole moment of the water molecule in the 
liquid phase of around 2.9 + 0.6 D (Badyal et al. 2000). The difference of around 1.1 + 0.6 D 
is thus an average induced dipole moment typical for liquid water. This approach works 
well for homogeneous surroundings, e.g. liquid water, but is less accurate for heterogeneous 
surroundings, e.g. interactions with ions or at surfaces, where the molecular induced dipole 
moment will be significantly different from the average value. 

Another example of effective force fields is united atom force fields, where functional groups 
are represented as pseudo-atoms. Most commonly, hydrogen atoms are suppressed in 
unpolar system, e.g. some polymers like polyethylene, where for example the methyl 
group -CH 3 is reduced to a one-centre Me group normally placed at the position of the C 
atom. Suppressing all hydrogen atoms will drastically reduce the number of interactions to 
calculate in eq. ( 3 . 3 . 59 ), and will speed up the calculation dramatically. 

Inorganic systems are in principle more complicated since they include heavier elements 
and d- and /-electrons would require higher order atomic multipoles as quadrupole and 
possibly octupole moments to describe the electrostatics. Many solid-state systems have, 
however, ionic bonds leading to that charge-charge interactions dominate. A simple 
example is the Born model for ionic solids, 


y _ y' _ did I + Aij 

1=1 7=7+1 477 Cpi?// R"j 

electrostatics repulsion 


(3.3.63) 


3.3.5 Summary of nonbonding terms 

The most important non-bonding energy terms in a force field for modeling intermolecular 
interactions as well as long-range interactions within molecules are 

• electrostatics 

• induction 

• dispersion 

• short-range repulsion 

In general for a force field, each energy term relies on a choice of a function form which 
in the best case are derived from quantum mechanics. Force fields also rely heavily on a 
set of atom-type parameters which needs to be parametrized from e.g. a set of atom-type 
parameters. So far we have given a more or less phenomenological introduction to force 
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fields. In section 3.4, the force field terms for intermolecular interactions are derived from 
quantum mechanical perturbation theory putting this section on a more firm basis. 

The advantage of force fields is that they present a rapid way of calculating interaction 
energies as well as the force acting on each particle in the system, and force fields are 
therefore used extensively in molecular dynamics and Monte Carlo simulations pf liquids 
where the interaction between many particles need to be calculated repeatedly. 


3.4 Intermolecular forces from quantum mechanics 


Perturbation theory may be adopted to link the force field approach to quantum mechan¬ 
ics (Buckingham 1967, Margenau and Kestner 1969, Stone 1996, Engkvist et al. 2000). For 
two molecules, A and B, where molecule A has n A electrons and Na nuclei and molecule 
B has tib electrons and Ng nuclei, respectively, the Hamiltonian may be divided into three 
contributions as 

je = je A +je B +VAB- 0 . 4 . 1 ) 

The molecular Hamiltonian of a molecule A is given from eq. (2.4.2) as 


tiA _ i n-A Na _ y 

4 /=!/=! ni 


i =1 


n A 

L • 

i—i, 

j=i +1 


+ 


’ij 


N a 

L 


ZiZj 


1=1, R IJ 
/=/+! 


(3.4.2) 
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and the part describing the interaction between the two molecules, Y AB > is given as 


ZiZj 


ft A N B _7 r Na ftB ft A ftB 1 Na N B 

^B = EE-^ + EE-^ + EEf+ EE * 

i=l /=1 r U 1=17=1 'Ij i=lj=l'ij 1=1 /=l n U 


(3.4.3) 


It is noted that we in eqs. (3.4.i)-(3.4.3) have assigned some nuclei and electrons to 
molecule A and the remaining nuclei and electrons to molecule B, and as compared to 
eq. (3.4.2) the only difference is that the summations have been restricted. Within the Born- 
Oppenheimer approximation the positions of the nuclei are known and it is in most cases 
trivial to assign a nucleus to a molecule. For the electrons, on the other hand, we have stated 
that ha electrons belong to molecule A and ns electrons belong to molecule B, which is an 
approximation. 

Regular Rayleigh-Schrodinger perturbation theory in section 2.8.1 is adopted on the inter¬ 
action between two molecules. Using multipole expansions (see section 2.9.1) and the long- 
range approximation for intermolecular interactions, it is demonstrated that second-order 
perturbation theory gives an electrostatic energy, an induction energy and the dispersion 
energy. 


3.4.1 The first-order energy 


Using the multipole expansion in eq. (2.9.35) leads to that the interaction operator in 
eq. (3.4.8) is given as 


OO 

Tab = Y 

m,n=0 


r-nm 

tl J? {m) j(m+n ) 

YYl\yi\ A,a\...a m AB,(X\...oc m+n B,oc m+ \...cc m+n 


(3.4.4) 


Adopting Rayleigh-Schrodinger perturbation theory, where Tab is regarded as a perturbation 
to the molecular Hamiltonians, ,76a and ,76%, in eq. (3.4.1), the first order correction to the 
energy becomes according to eq. (2.8.9), 

£ 0 ] = ’ (3.4.5) 


where the ground state is considered. The molecular wavefunctions are supposed to 
be known, e.g. J&aWaj) = £ A,iWA,i), and the molecular wavefunction if/ a, i is correctly 
antisymmetrized with respect to the ha electrons of molecule A. The wavefunction of two 
molecules may be written in terms of the molecular wave functions as a generalized Heitler- 
London wavefunction (Margenau and Kestner 1969), 

I^[, 0) >=^4bI^a ) o'I/'b,o), (3.4.6) 

where s7ab is an operator that ensures that the total wavefunction is correctly antisym¬ 
metrized. In the long-range approximation, it is assumed, however, that the zeroth-order 
wavefunction is given as a direct product of the two molecular wavefunctions, 

|t/4 0) > « \ysA,oV/B,o) • (3.4.7) 

The operator, s7ab> will construct a set of correctly anti-symmetrized orbitals for the entire 
complex. This effect will only be substantial if the wavefunctions for the separate molecules 
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overlap, Sab = (Vsa,o\Vsb,o)> which declines exponentially with the distance between the two 
molecules. Since y/A ,o is a function of only the electronic coordinates of molecule A, and 
likewise t (/b,o is a function of only the electronic coordinates of molecule B, the first-order 
energy becomes by using eq. (3.4.4), 


£ 


(i) 

o 


+ 


oo 


(- 1 ) 
mini 

OO 


E 

m,n =0 


m 




( m+n ) 


Hn) 


Mm+l’-ttm+n 


Wb,q) 
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(m) 


Am+ri) 


M 


in) 


m,n =0 


m \ n \ A } oc\...0L m AB f (x\...oc m + n B,oc m+ \...oc m + n 


qATf B q B + q A T^ ] Ba id B ,a + a/ 3 Q B , aj6 +... 

qA,a ^46,a qB - RAa ^AB,al3^ B ’P ~ 2^ A,a ^AB,a/3r^ B >^T + ‘" 

2^ A ’ a P^AB,al5^ B + 2^ A,a ^A B,a^Y^ B A + ^^ A ’ a P^AB,a/3jS^ B ’T 5 + ■■■ > (3.4.8) 


where a molecular moment, , is defined in eq. (2.9.50). The electrostatic energy in 

eq. (3.4.8) may be rewritten in terms of the quadrupole moment in eq. (2.9.32)) as 


,(i) 


qA. T AB 4 B + C I A ^AB,a^ B ’ a + 3 C ! a ^AB,ap^ B > a P + '" 

T-(l) 2) 1 ^(3) a 

^Aa 1 AB,a°l B ^Aa i A B,a^ B ’P 3 ^Aa 1 A B,apy UB ’Pr + ■■■ 


+ ^ A ’ a P T AB,ap £ f B + ^ A ’ a P T AB,ap r ^ B ’ r+ ^ A ’ a P T AB,ap r S^ B A s + --- ’ ( 3 - 4 - 9 ) 

In figure 2.6, it is noted that the molecular electric moments are calculated around the local 
origin of each molecule. Adopting the definition of the electrostatic potential in eq. (2.9.1), 
the electrostatic potential at molecule A becomes 


q> A - T ab VB + T AB,a^ B ’ a + 2 T AB,ap ^ B ^P + ■ ■ ■ 

CO | 

_ + — jin) j^jin) 

Yi | AB } 0C\...0Cyi B } 0C\...0Cyi 
n— 0 n • 


(3.4.10) 


and the electrostatic potential at molecule B becomes, 


<Pb 


<7a T ab Ha, a Tab, a + ~QA,ap T ABaB + 


-d) 


( 2 ) 


CO 


E 

n =0 


_i n 

— M (n) 

jq] A y OC\...OC n 


r rin) 

AB yOc\...oc n 


(3.4.11) 


where the difference in sign for the dipole moment terms again arises from the definition 
of the distance vector Rab in figure 2.6. The mth derivative of the electrostatic potential at 
molecule A with respect to a displacement is given as 


im) 

VA,a\. 


oo 1 

= y -t 

n ! 1 

n =0 m 


{m+n) 


M' 


in) 


AByGC\...OC m+n B,cc m + \...oc m + n 


(3.4.12) 
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and correspondingly 


(m) 



_i n 

j{m+n ) 

jq] A,OCm+\'"Clm+n 


(3.4.13) 


The first-order energy may thus be written as 


oo _ i m 

e w = y —-—M (m ^ (o {m) 

0 yyi\ A,OC\...OL m ' A,OC\...OC m 

m =0 ,n ' 

OO 1 

= y —w (m] M {m] 

rn\'B t OC\...OCm B t OC\...OC m 
m =0 " l - 


(3.4.14) 


The electrostatic energy may thus be calculated either as the interaction between the electric 
moments at the molecules or as the interaction between the electric moment and the Taylor 
expansion of the electrostatic potential at one of the molecules. 


3.4.2 The second-order energy 

For the interaction between two molecules in the long-range limit, where it is again 
assumed that the zeroth-order wave function is given as the direct product of two molecular 
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wavefunctions as in eq. (3.4.7), the second-order energy in eq. (2.8.18) becomes 

( 2 ) ^ (VA,0V'B,0\yAB\y'A,uV'B,v)(V'A,uV'B,v\yAB\VA,0VB,0) „ , , „ 

£ 0 = y - (3.4.15) 

u+v— i £a,o-£a,u + £b,o~£b,v 

where the notation | y/^) = |i/m,oV / b,o)> £q 0) = £a,o + £b ,o and Ep ] = Ea,u + £b,v have been 
adopted. Using the multipole expansion of Tab in eq. i (3.4.4), the second-order energy 
becomes, 
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(3.4.16) 


This energy is divided into three terms that are treated separately: a) u = 0, v = 1,... ,oo, b) 
u = l,...,oo, v = 0,c) u = l,...,oo, v = l,...,oo. The first term {u = 0, p = l,...,oo) becomes 
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(3.4.17) 


which is analogous to the expression obtained in eq. (2.9.60) for a molecule in an external 
electrostatic potential. The leading term {n = 1, q = 1) becomes 


£ 0,a ~ ~2 aB ’ a P EB ’P EB ’ a 


(3.4.18) 


which as in eq. (2.9.63) is identified as a induction energy. The second term, {u = l,...,oo, 
v = 0), becomes 
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If eq. (3.4.17) gives the energy of the polarization of particle B, eq. (3.4.19) gives the 
contribution to the polarization of particle A. Again the leading term (m = 1, p = 1) is given 
as 


£ 0,b~ 2 

and the total induction energy, £^ a+b > becomes 


(3.4.20) 
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(3.4.21) 


The result obtained in eq. (3.4.21) is not consistent with the classical expression for the 
induction energy given in eq. (3.4.24) since the contribution to the electric field from the 
induced dipole moments is missing in eq. (3.4.21). It has been demonstrated that the 
additional term in eq. (3.4.24) may be obtained from higher order energy terms in the 
perturbation expansion (Stone 1989). 

The third term (u = 1,..., oo, v = l,...,oo), which will be identified as the dispersion energy, 
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(3.4.22) 


is not as straightforward as the induction energy to separate into term that are dependent on 
the properties of the separate molecules. The nominator is readily separated into a product 
of molecular properties, but in order to resolve the denominator, the following identity is 
used, 
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The dispersion energy is thus rewritten as 
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where the notation coj <n = e„ - £/ o has been introduced. The frequency-dependent dipole- 
dipole polarizability, a a p(to), is defined in eq. (2.9.75) leading to 
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If again, only the leading term in the second order energy, (m = 1, n = 1, p = 1, q = 1), is 
included, the dispersion energy is given as 
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(3.4.26) 


If the frequency-dependent polarizability is approximated according to the Unsold approx¬ 
imation (Unsold 1927) 
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where w is a kind of average excitation energy of the molecule, the dispersion energy 
becomes 



1 

4 


(OaWb 
to a + 



(3.4.28) 


Finally, if isotropic polarizabilities, a a p = a Sap, are regarded, the expression for the 
dispersion energy is reduced to the well-known London formula, 



Ce^AB 

d6 

K AB 


(3.4.29) 


where C^,ab is an atom-pair parameter depending on the polarizabilities. 

To conclude, the second order energy provides both the induction and dispersion energies. 
Both of these energies describes a response to the presence of the other molecule. Since a 
system always responds such as the energy is minimized, both the induction and dispersion 
energies give attractive contributions to the interaction energy. 
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Exercises 


Ex. 3.1 — Derive the result for the two-particle system of the EEM in eq. (3.3.14) and from its 
general equations in eq. (3.3.10) and (3.3.11), 


—< 7 / = < 7 / 


tj-ti 

r ll + i 1 j-2Tf ) J > 


Also explain why this is a reasonable model for 

i) having one set of atom-type parameters for each element describing the atomic charges 
in both e.g. CO and O 2 . 

ii) describing the geometry-dependence of atomic charges. 

Finally, what is the fundamental flaw in the EEM (seen from inspection of the two-particle 
equation at R — 00 )? 


Ex. 3.2 — With a self-energy, T { ff = rij, the Coulomb interaction energy for a set of N atoms 
may be written as 

z I,K 

including also the I = K term. We change to charge-transfer variables, 

N N 

= qij and q k = Y^ qKL 

J L 

where we imply that qn = 0 (i.e. there is no charge transfer to "itself") and that charge 
conservation is imposed by qu = -qji. Show that V can be rewritten in terms of the charge- 
transfer variables qij in a "symmetrized" way (i.e. with respect to the subscripts I and J as 
well as K and L, respectively) so that 


N 


V: 


~ L 

8 I,hc,L 


qu( 


r r{ 0) . y 

1 IK “ r 1 JL 


( 0 ) 


-T 


( 0 ) 


IL 



Ex. 3.3— Consider the interaction energy between methane and the fluorine ion, F _ . In 
terms of a one-centre multipole expansion for each molecule (ion), which is the leading 
term? (Hint: model methane with atomic charges and use carbon as the local origin in 
methane for the calculation of its electric moments) 

Ex. 3.4 — Consider a system with two atomic charges with magnitude + q placed in (0,0, ±z ) 
and an isotropic point polarizability a placed in (0,0,0). What is the interaction energy 
(electrostatics and polarization)? Is the interaction energy pair-wise additive? 
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electronegativity equalization model, 85 

electrostatic potential, 31 

elementary charge, 7 

empirical force field, 74 

energy constraint, 79 

energy restraint, 80 

exchange integral, 18 

exchange operator, 45 

exchange-correlation functional, 61 

exchange-repulsion energy, 95 

expectation value, 9 

Fermi's golden rule, 29 
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